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1  Introduction  and  Scope  of  the  Problem 


1 . 1  Foliage  above  a  Rough  Surface:  An  Approach  Rationale 

When  studying  scattering  from  the  combination  of  a  foliage  layer  above  a  rough  surface, 
one  must  deal  with  the  scattering  from  the  foliage,  the  rough  surface,  and  the  interaction 
between  the  two.  Of  these  three  components  to  the  fundamental  scattering  problem,  the 
foliage  scattering  and  the  interaction  scattering  are  definitely  the  most  difficult  to  analyze. 

Foliage  by  itself  presents  a  real  challenge  to  the  analyst  because  of  the  many-body  aspects 
of  the  problem.  In  addition,  modeling  foliage  scattering  from  first  principles  is  further 
complicated  by  the  fact  that  the  scatterers  are  irregular  at  best  and  generally  ill  defined.  From 
an  electromagnetic  point  of  view,  leaves  do  not  all  look  to  be  identical  and  twigs,  branches, 
and  limbs  conform  to  no  particular  shape!  Thus,  except  for  very  low  and  very  high  frequency 
limits,  it  is  almost  impossible  to  compute  the  scattering  pattern  of  the  basic  "constituents"  of 
foliage.  We  do  know  a  bit  about  the  typical  volume  fraction  of  foliage  and  how  this  is 
partitioned  between  leaves  and  the  woody  components,  and  we  have  some  idea  of  the  range 
of  complex  dielectric  constant  variation  for  wood  and  leaf  materials  [1],  Yet  another 
unknown  is  the  variation  of  foliage  density  with  depth  into  a  canopy.  Finally,  even  though 
foliage  is  quite  frequently  classified  as  on  the  edge  of  being  a  volumetrically  sparse  medium, 
this  does  not  mean  that  there  is  a  lack  of  strong  electromagnetic  interactions  between  the 
various  scattering  components,  e.g.,  twigs,  branches,  leaves,  etc.  Furthermore,  for  European 
forests  that  have  been  well  managed  and  not  logged  (also  called  old  growth  forests),  the 
volume  fractional  density  of  the  biomaterial  may  be  as  large  as  5%. 

When  dealing  with  independently  scattering  objects  and  scattering  from  rough  surfaces,  it 
is  possible  to  convert  single-frequency  models  of  the  individual  scattering  cross  sections  (for 
the  discrete  objects)  and  the  scattering  cross  section  per  unit  area  (for  the  extended  surface 
scattering)  into  models  for  the  incoherent  time-dependent  scattered  waveform  produced 
under  pulse  illumination.  For  strongly  interacting  individual  scatterers,  this  simplification  is 
not  usually  possible,  the  reason  being  that  it  is  not  sufficient  to  know  that  the  scatterer  is  in  a 
volume  because  its  location  within  the  volume  must  also  be  known. 
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What  can  be  done  then  to  resolve  this  dilemma?  First,  we  know  that  we  can  always  extract 
an  effective  “scattering  cross  section  per  unit  illuminated  volume,  ov”  from  airborne  radar 
data,  e.g.,  the  scattered  power  Pr  received  by  a  pulsed  radar  is  given  by 


Pr 


1  (4tt)2R4 


(1) 


where  the  effective  scattering  cross  is  given  by 

c  =  (cT/2)(R®a2)(R0tlK  (2) 


Solving  (1)  and  (2)  for  av  yields 

„  .£ _ <4*>2r2 _  (3) 

P,  G2(0i,*i)(cT/2)(4>!l2eel) 

In  the  above  equations  Pr  is  the  received  or  scattered  power,  Pt  is  the  transmitted  power, 
is  the  two-way  antenna  gain  in  the  indicated  direction,  c  is  the  speed  of  light,  T  is 
the  pulse  length,  R  is  the  range  to  the  volume,  Oaz  is  the  antenna’s  azimuthal  beamwidth,  and 
©ei  is  its  elevation  beamwidth.  Of  course,  some  of  these  factors  are  a  function  of  time 
indicating  what  portion  of  the  volume  scatterers  are  being  illuminated  by  the  incident  pulse 
waveform  as  it  passes  through  the  scattering  medium.  Equation  (3)  is  actually  an 
approximation  in  that  one  should  integrate  over  the  volume  bounded  by  the  antenna  gain 
pattern  weighting  and  the  pulse  width  extent.  The  important  point  is  that  incoherent  power 
waveform  data  from  a  pulsed  radar  can  be  converted  into  an  effective  "scattering  cross 
section  per  unit  illuminated  volume".  It  should  be  noted  that,  within  the  resolution  limits 
imposed  by  the  radar  pulse  width  and  antenna  beamwidth,  the  av  extracted  from  the  data  may 
be  a  function  of  the  slant  path  distance  (R)  into  the  medium. 

The  next  step  in  the  process  is  to  compare  these  measured  data  with  our  models  and  it  is 
here  that  things  become  difficult.  First,  as  noted  above  there  are  no  tractable  "first- 
principles"  models  for  propagation  and  scattering  by  a  foliated  environment  under 
conditions  of  strong  interaction  (multiple  scattering).  Secondly,  those  models  that  claim  to  be 
"more"  applicable  to  such  an  environment  usually  are  deeply  imbedded  with  involved  and 
tedious  computations  whose  physical  meaning  is  marginal  at  best.  This  quandary  means  that 
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a  decision  must  be  made  as  to  how  to  best  proceed,  i.e.,  an  approach  rationale  must  be 
developed. 

The  approach  rationale  followed  in  developing  the  model  presented  in  this  report  is  as 
follows.  First,  it  is  well  known  that  the  end  users  of  such  models  are  extremely  skeptical  of 
any  model  that  does  not  contain  some  measurements  in  its  development.  That  is,  they  are 
concerned  that  the  model  be  designed  so  that  it  is  capable  of  at  least  reproducing  known 
measurements.  Consequently,  we  felt  it  essential  to  involve  measured  data  in  our  model.  The 
second  element  of  this  model  is  based  upon  the  realization  that  it  is  possible  to  develop  a 
somewhat  general  model  for  media  that  is  not  too  strongly  interacting  and  this  model  may  be 
"matched"  to  data  to  determine  the  actual  parameters  that  are  embedded  in  the  model  and, 
perhaps,  to  extend  it  beyond  its  known  range  of  validity.  In  short,  the  model  parameters  can 
be  determined  by  matching  actual  measured  data  to  the  scattering  results  predicted  by  the 
model.  It  should  be  noted  that  such  an  approach  avoids  long  and  questionable  computations 
based  on  one's  estimate  of  what  actually  causes  the  scattering  and  how  it  does  this.  The 
reason  for  avoiding  such  computations  is  very  simple  -  there  is  no  way  to  estimate  how 
applicable  they  will  be  since  the  accuracy  of  the  overall  model  is  unknown.  By  matching  the 
model  to  data,  we  are  in  effect  extending  the  model's  accuracy  through  the  use  of  "effective 
parameters"  that  are  "calibrated"  by  the  data.  In  summary  then,  our  approach  has  been  to  do 
the  best  we  can  to  develop  a  model  that  is  accurate  but  not  overly  (computationally)  detailed, 
match  it  to  measured  foliage  scattering  data  to  generate  the  model  parameters  ("effective 
parameters"),  and  then  investigate  the  accuracy  of  extending  the  model  to  other  situations 
using  the  existing  "effective  parameters". 

The  remainder  of  this  report  details  efforts  to  develop  the  "best"  model  that  we  can  for  the 
foliage  and  its  interaction  with  the  terrain  using  the  above  rationale.  We  emphasize  that  while 
it  may  not  be  esthetically  pleasing  to  have  to  resort  to  measurements  to  truly  complete  the 
model,  this  guarantees  that  the  model  parameters  derived  from  the  data  will  be  valid  and,  in 
fact,  may  compensate  for  certain  deficiencies  in  the  model.  Given  the  end-goal  of  predicting 
foliage  scattering  and  penetration  and  trying  to  come  up  with  new  methods/techniques  to 
penetrate  foliage,  such  an  approach  as  this  may  be  the  most  appropriate! 
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1.2  The  Major  Contributors  to  the  Return  Waveform 


Ideally,  in  estimating  the  returned  signal  from  foliage  covered  terrain,  a  scattering  model 
will  include  self  and  mutual  interactions  among  the  constituent  components.  In  the  following 
sections,  a  brief  description  of  each  significant  interaction  will  be  presented.  The  presentation 
in  these  sections  will  follow  the  integral  equation  approach. 

We  begin  by  considering  the  incident  field,  Emc ,  that  exists  in  free  space  in  the  absence 
of  the  foliage  and  the  surface,  see  Figure  1.  With  introduction  of  a  scatterer,  i.e.  foliage,  the 
total  field  is  found  to  be  a  superposition  of  the  incident  field  E,nc  and  the  field  scattered  by 
the  scatterer,  Es. 

Etotal  =  Einc  +  Es 

Consequently,  our  first  task  will  be  the  construction  of  the  scattered  field  produced  for  the  ith 
scatterer,  E- ,  of  the  N  objects  which  comprise  the  volume  of  scatterers. 


Figure  1:  The  Incident  Field 


4 


1.2.1  Scattering  from  a  Volume  of  Discrete,  Closed  Body  Scatterers 

Each  component  of  the  foliage  (leaves,  twigs,  branches,  etc.)  will  scatter  energy  from  the 
incident  field.  Isolated,  each  scatterer’s  effect  can  be  assessed  using  an  integral  equation 
approach  that  leads  to  a  method  of  Moments  (MOM)  formulation.  Consequently,  the 
scattering  solution  is  found  in  terms  of  a  vector  scattering  pattern  with  respect  to  the  incident 
field  and  the  incident  field’s  angle  of  arrival.  This  scattering  pattern  is  then  used  to  construct 
the  scattered  field  due  to  an  object.  Typically,  the  scattering  pattern  is  derived  using  the  far- 
field  approximation  and  assuming  an  incident  plane  wave. 

Considering  a  volume  of  N  scatterers  without  mutual  interaction,  the  scattering  patterns 
of  the  individual  members  can  be  used  to  construct  the  volume  response  in  the  direction  of 
observation.  When  the  location,  shape  or  size  of  the  individual  scatterers  are  uncorrelated,  a 
simple  summation  of  the  scattered  power  waveform  from  each  scatterer  is  possible.  This 
describes  the  results  of  the  single  scatter  theory  and  first  order  multiple  scattering,  see 
Appendix  A.  First  order  multiple  scattering  plays  a  central  role  in  the  model  proposed  in  this 
report.  It  is  presented  in  a  general  form  in  Appendix  A  and  a  restricted  form  in  Chapter  2. 

If  multiple  scattering  within  the  volume  is  expected  to  play  an  important  role,  interactions 
among  the  discrete  scatterers  must  be  considered.  An  exact  solution  would  consider  not  only 
the  scattering  from  each  individual  object,  but  the  full  interaction  between  them.  One  method 
of  solution  that  accounts  for  full  interaction  between  N  objects  is  the  solution  of  the 
associated  N-coupled  integral  equations.  Using  the  equivalence  principle,  the  scatterers  may 
be  replaced  by  equivalent  currents  that  radiate  in  free  space.  In  addition  to  the  requirement 
for  more  computing  power  than  is  commonly  available,  the  exact  solution  for  the  N-coupled 
integral  equations  for  propagation  through  the  foliage  would  require  an  abundance  of  detailed 
data  and  a  number  of  realizations  to  create  acceptable  averages.  Consequently,  some 
approximations  must  be  made. 

The  reduction  of  the  problem  for  backscatter  and  propagation  through  the  random  foliage 
can  be  achieved  in  several  ways,  depending  on  the  level  of  interactions  to  be  included.  For 
propagation  of  the  mean  field,  approximations  in  the  literature  include,  but  are  not  limited  to 
the  following 
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•  Single  scattering  theory  (see  Appendix  A) 

•  First  order  multiple  scattering  theory  (see  Appendix  A) 

•  The  Foldy-Lax-Twersky  Integral  Equation 

Higher  order  moments  of  the  field  can  also  be  handled  by  forming  the  proper  moments  in 
the  first  two  cases.  These  higher  order  moments  are  necessary  in  order  to  account  for  the 
propagation  of  pulses.  Twersky,  among  others,  has  developed  integral  equations  that  describe 
the  higher  order  moments  which  include  some  level  of  multiple  scattering  [Ishimaru,  1997]; 
various  approximate  solutions  exist  for  these  equations.  In  addition,  there  are  also  hybrid 
techniques,  such  as  the  Distorted  Wave  Bom  Approximation  (DWBA)  [Lang,  1981],  Lang 
used  the  Foldy-Lax-Twersky  Integral  Equation  in  order  to  establish  the  mean  field;  he  then 
used  the  single  scattering  theory  for  particles  immersed  in  an  equivalent  media  derived  from 
the  mean  field  to  find  the  second  moment.  Our  approach  is  similar  to  this  DWBA  approach. 

Whether  single  scatter  theory  or  multiple  scatter  theory  is  used,  we  will  construct  a 
composite  scattered  field,  E, ,  due  to  the  volume  of  scatterers,  see  Figure  2.  In  the  integral 
equation  formulation,  the  scattered  field  results  from  an  induced  current,  Jsn,  on  each  of  the  N 
scatterers,  see  Figure  2. 


Figure  2:  The  total  incident  field  with  respect  to  the  surface 

Initially,  we  present  a  derivation  of  the  scattered  field  from  the  volume  using  a  reduced  form 
of  the  radiative  transfer  approach  that  results  in  a  form  of  the  first  order  multiple  scattering 
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result;  see  Chapter  4  and  Appendix  A.  The  limitations  of  this  result  are  discussed  at  the  end 
of  Appendix  A  in  the  context  of  the  first  order  multiple  scattering  results. 

1.2.2  Scattering  from  a  Rough  Surface 

To  model  the  multiple  scattering  that  takes  place  along  the  surface,  there  are  several 
techniques  available.  Numerical  implementations  typically  use  the  integral  equation 
formulation.  Given  a  field  incident  on  a  statistically  rough  terrain,  the  integral  equation 
technique,  typically  numerically  implemented  with  the  Method  of  Moments  (MOM),  can 
yield  exact  results  for  a  given  realization.  Average  results  for  a  collection  of  realizations  are 
found  using  Monte  Carlo  methods.  A  solution  method  for  the  MOM  formulation  that  is  of 
interest  in  this  report  is  the  Method  of  Ordered  Multiple  Interactions  (MOMI).  It  will  be 
explored  in  Chapter  4  and  is  used  to  verify  some  assumptions  in  the  model  developed  in  this 
report.  Analytical  results,  including  Kirchhoff  and  perturbation  approximations,  may  be 
useful  in  an  analytically  reduced  integral  equation  formulation  discussed  in  Chapter  4. 

A  time  dependent  analytical  approach  for  the  calculation  of  the  incoherent  power 
waveform  from  an  extended  rough  surface  that  is  consistent  with  the  single  scatter  approach 
is  the  Impulse  Response  Method.  This  method  is  derived  under  the  assumption  that  there 
exist  a  continuum  of  scattering  facets  on  the  surface  that  reflect  a  radar  power  waveform 
[Brown,  1977].  Under  certain  assumptions,  the  return  power  from  each  properly  oriented 
surface  facet  enters  into  a  summation.  The  number  of  these  properly  oriented  facets  per  unit 
area  of  the  surface  defines  a  cross  section  per  unit  area.  Chapter  2  presents  a  brief  discussion 
of  the  impulse  response  method  for  calculating  the  scattering  from  terrain  in  free  space,  i.e. 
no  foliage  cover. 
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1.2.3  Interaction  between  the  Foliage  and  the  Rough  Surface 

Regardless  of  the  modeling  method  used  for  the  terrain  and  the  volume  scattering 
individually,  there  will  be  an  additional  source  of  multiple  interactions  (or  multiple 
scattering):  the  interaction  between  the  volume  and  the  surface.  Once  the  field  scattered  from 
the  surface,  ,  due  to  the  currents  induced  on  the  surface,  Js ,  is  calculated,  it  can  act  as  an 
additional  field  incident  upon  the  foliage  on  its  path  back  through  the  foliage  to  the  radar.  See 
the  field  in  Figure  3. 


Js 

Figure  3:  Surface  to  Foliage  Interaction 

This  single  passage  (from  foliage  to  surface,  and  back  to  foliage)  does  not  account  for 
the  full  interaction  between  the  current  induced  in  the  foliage  and  the  current  on  the  surface. 
This  single  passage  approximation  to  the  interaction  between  the  foliage  and  the  surface 
represents  a  single  interaction:  the  foliage-scattered  field  that  creates  the  surface  currents  is 
due  only  to  the  field  incident  from  the  radar.  This  approximation  is  explored  in  Chapter  4 
through  a  comparison  with  the  exact  results  for  a  single  scatterer  above  a  rough  surface.  A 
full  interaction  formulation  requires  that  the  currents  on  the  surface  and  on  each  scatterer  in 
the  foliage  be  coupled.  Additional  interaction  terms  would  include  corrections  to  the  foliage 
currents  due  to  the  surface  scattered  field  that,  in  turn,  will  produce  corrections  to  the  surface 
currents:  an  infinite  series  of  these  corrections  will  produce  the  full  interaction  results. 
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Alternatively,  a  coupled  integral  equation  formulation  relating  the  induced  currents  will  also 
produce  the  full  interaction  result;  this  method  for  verification  is  explored  in  Chapter  4. 

Assuming  that  the  integral  equation  approach  is  followed,  the  solution  for  the  passage 
from  the  surface  through  the  foliage  can  be  formulated  using  the  equivalence  principle.  The 
surface  current,  Js,  is  permitted  to  radiate  in  free  space  and  the  resulting  field  acts  as  an 
additional  incident  field  with  respect  to  the  foliage.  Hence,  the  surface  scattered  field.  Eg , 
will  induce  a  corrective  current,  Jssn ,  on  the  nth  scatterer  which  must  be  vectorially  added  to 
the  previous  current.  This  current  radiates  a  second  field  scattered  from  the  foliage,  Egf ,  in 

addition  to  the  scattered  field  due  only  to  the  foliage,  E* .  See  Figure  4.  This  is  a  first 
approximation  to  the  foliage-surface-foliage  interaction. 

A  second  order  correction  to  the  surface  current  will  treat  the  incident  field  on  the 
surface  as 

ginc  on  surface  _  ginc  gs  gs 

Consequently,  a  new  surface  current  Js  is  found.  This  current  will  produce  a  new  value  for 
the  surface  scattered  field,  Eg ,  and  a  new  value  for  the  field  incident  to  the  foliage  from  the 
surface:  Continuing  this  process  of  iteration  will  produce  the  full  interaction  result. 

Alternatively,  like  the  first  passage  through  the  vegetation,  this  second  field  scattered 
from  the  foliage,  Esfef ,  can  be  found  using  the  single  scatter  approximation.  This  result  may 

also  be  iterated,  correcting  the  surface  currents  producing  scattered  fields.  However,  this 
result  will  suffer  the  limitations  of  the  single  scatter  theory  as  discussed  in  Appendix  A. 
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Figure  4:  the  second  order  approximate  scattered  field  from  the  foliage  and  surface 

combination 


1 .3  Goals  of  this  Research 

The  goals  of  this  research  are  to  produce  a  numerically  efficient  code  which  can 
incorporate  measured  data  for  calibration  and  accurately  reproduce  the  general  trends  of  an 
average  returned  waveform  from  terrain  and  foliage  with  similar  statistics  and  constituents  as 
the  calibration  data.  Numerical  efficiency  is  best  served  through  the  use  of  the  impulse 
response  approach,  which  casts  the  returned  waveform  into  a  series  of  convolutions  and  uses 
empirically  derived  parameters.  Through  its  relation  to  first  order  multiple  scattering  theory, 
this  model  was  found  to  incorporate  some  limiting  assumptions.  These  assumptions  are 
discussed  in  Appendix  A.  In  order  to  test  these  assumptions,  the  generalized  form  of  the 
radiative  transfer  result,  first  order  multiple  scattering  model  has  been  expanded  to  simulate  a 
single  scatterer  over  a  rough  surface.  Some  of  these  assumptions  are  quantified  to  a  certain 
extent  in  the  later  sections  of  Chapter  4  by  comparison  with  the  “exact”  numerical  results. 
Although  the  simplest  model,  the  radiative  transfer  result,  described  in  Chapter  3,  may  not  be 
totally  adequate,  at  least  some  aspects  of  the  convolutional  formulation  may  be  validated.  In 
addition,  the  convolutional  form  has  been  retained  through  the  use  of  the  more  general  first 
order  multiple  scattering  theory  found  in  Chapter  4  and  Appendix  A.  This  development 
represents  the  major  thrust  of  the  work  remaining  to  be  done  in  order  to  verify  the  fully 
convolutional  result. 
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In  addition  to  verification  via  the  first  order  multiple  scattering,  alternate  formulations 
that  incorporate  more  multiple  scattering  elements  are  discussed  in  the  later  sections  of 
Chapter  4,  including  a  reduced  integral  equation  result.  Finally,  it  is  postulated  that  the 
results  generated  by  a  first  order  multiple  scattering  result,  using  Foldy-Lax-Twersky  integral 
equation  for  the  mean  field  to  find  the  effective  media,  will  adequately  approximate  the 
results  of  Lang  [1981].  Consequently,  a  comparison  with  Lang’s  results  from  the  Distorted 
Wave  Bom  Approximation  (DWBA)  will  be  necessary.  This  effort  is  important  to  better 
understand  the  limitations  of  the  model  and  for  model  verification. 


ll 


2  Estimating  the  Isolated  Surface  Return  Power 

2.1  Integral  Equation  Formulation  of  Rough  Surface  Scattering 
and  the  Method  of  Ordered  Multiple  Interactions 


The  integral  equation  governing  both  the  TE  and  TM  polarizations  in  the  2-D  scalar 
problem  has  been  derived  in  many  sources  including  [Ishimaru,  1994],  by  many  different 
techniques,  such  as  equivalence  and  the  use  of  Green’s  Identities.  Green’s  second  identity 
is  given  by 


f(r)-f'(f)+j 


dn' 


dn' 


r  €  V 


(1) 


In  (1),  V  represents  a  certain  volume  in  space  surrounded  by  the  closed  surface  S  and  S* 
as  shown  in  Figure  5. 


Figure  5:  Problem  Geometry  for  the  Derivation  of  Boundary  Integra]  Equations 

Referring  to  Ishimaru  [1994],  an  intermediate  result  derived  by  these  methods  is  our 
starting  point;  we  start  with  the  following  integral  equation  relating  the  total  scalar  field, 
f(r) ,  at  the  observation  point,  r ,  to  the  incident  field,  f '(r ) : 


f(f)  =  2fi(f)+J 


an' 


an' 


dl',  reC 


(2) 
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where  the  contour  integration  is  taken  over  a  boundary  enclosing  the  source  and 
observation  points  consists  of  a  hemisphere  at  infinity  and  a  contour  along  the  rough 
surface;  see  Figure  5. 

The  Green’s  function,  G(r,r'),  and  its  normal  derivative  are  chosen  such  that  their 
contribution  is  zero  at  infinity;  hence,  only  the  integral  over  the  surface  remains.  This 
integration  over  the  surface  is  reduced  in  its  support  by  limiting  the  illuminated  region  to 
only  a  portion  of  the  surface  by  use  of  a  tapered  beam.  A  more  detailed  description  of  this 
process  can  be  found  in  many  sources  including  Kapp  and  Brown  [1996]. 

For  1-D  surfaces,  the  contour  length  may  be  projected  onto  the  x-axis.  This  reduces 
the  integration  to  an  integral  over  one  Cartesian  coordinate.  Hence,  employing  the 
transformation 

dr  =  Vl-Kx(0  dx'  (2) 

we  can  construct  the  governing  integral  equations  in  rough  surface  scattering.  First,  the 
Electric  Field  Integral  equation  (EFIE),  can  be  derived  directly  from  (1)  by  enforcing  the 
following  boundary  condition  on  the  perfectly  conducting  surface:  f(r)  =  Ey(r)=0.  This 

results  in  the  following  first  kind  integral  equation  applicable  for  TE  polarization. 

Ei <f)  =  G(f , nji +£(*')  dx-  (3) 

—oo 

In  deriving  the  Magnetic  Field  Integral  Equation  (MFIE),  the  following  boundary 
condition  is  enforced:  the  normal  derivative  of  the  tangential  magnetic 
field, 0f(F)/3n’  =  3Hy(f')/3n’ ,  is  zero  on  the  surface.  This  results  in  the  following 

second  kind  integral  equation  applicable  for  TM  polarization. 

H,(f )  =  2H;(?)  +  2  +  CiiMd*'  (4) 

—CO 

In  order  to  express  the  equation  governing  the  TE  polarization  in  the  form  similar  to  the 
MFIE,  we  take  the  normal  derivative  of  both  sides  of  (3)  along  the  unit  normal  n  defined 
at  the  observation  point  r .  Then,  we  eliminate  the  weak  singularity  of  the  normal 
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derivative  of  the  Green’s  function  through  a  limiting  process  [Ishimaru,  1994].  This 
yields  the  following  second  kind  integral  equation  for  the  TE  polarization 


dEy(r)  dE'y(f)  %3Ey(r)  3G(rtr’) 

J  Pm' 


/l  +  Cx(x')dx' 


The  discretized  versions  of  the  above  equations,  when  properly  sampled,  yield  large,  full 
matrices  that  scale  as  the  number  of  unknowns  squared.  Scattering  from  a  rough  terrain, 
formulated  with  this  integral  equation  approach  typically  was  limited  to  small  surfaces  or 
narrow  incident  beams  due  to  the  matrix  storage  and  inversion  requirements  of  the 
conventional  method  of  moments  (MOM).  Solving  the  integral  equations  numerically  via 
the  Method  of  Ordered  Multiple  Interactions  (MOMI),  however,  has  reduced  this 
computation  time  and  storage  without  approximation  [Kapp  and  Brown,  1996].  Rewriting 
the  above  form  of  the  second  kind  integral  equations  as 


J(x)  =  J!(x)+  jK(x,x')J(x')dx' 


where  J(x)  is  the  unknown  surface  current,  K(x,x')  is  the  kernel  or  the  propagator,  and 
J‘(x)  is  known,  the  Kirchhoff  current.  Although  the  domain  of  integration  D  is  infinite 
by  design,  it  can  be  made  finite  with  the  use  of  the  appropriate  tapered  incident  field.  For 
the  TE  and  TM  cases,  respectively 

6E‘(x,z)|  .  , 

J'OO  =  2 ■  -gn—  U,,) .  J'(x)  =  2H',(x, z)|,.{w 

J(x')  =  ,  J(x')  =  Hy(x',z')| 

=  K(x.x’)  =  2 8G^’,X:i V1  +  £ OO 

dn  On 

After  discretizing  the  resulting  second  kind  equation  and  expressing  it  in  a  vector-matrix 


J  =  J1  +  PJ 
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In  (7)  both  J  (unknown)  and  J'  (known)  are  vectors  and  P  is  a  square  propagator  matrix. 
The  discretization  is  commonly  carried  out  by  taking  values  of  surface  height,  current  and 
propagator  at  the  uniform  grid  {xm}  of  N  discrete  points  separated  by  the  spacing  Ax .  In 
this  case,  the  mth  element  of  each  one  of  the  above  vectors  and  the  (m,n)th  element  of  the 
propagator  matrix  are  given  by 

Jm=J(xm)>  JL=Ji(xm)  and  P™  =P(xm,xn)Ax. 


The  off  diagonal  elements,  Pmn  with  imm,  of  the  discretized  propagator  matrix  P  are 
given  by 


dn 


dn' 


where  xm  -  (2m  -  l)Ax  -  N  Ax/2,  m  =  1,..., N  (observation  point  on  the  surface) 
xn  =  (2 n  -  l)Ax  -N  Ax/2 ,  n  =  1  (source  point  on  the  surface) 


for  the  TE  case  and  the  TM  case,  respectively.  The  diagonal  elements  (usually  called 
“self  terms”),  however,  require  special  treatment  and  are  given  by,  [Toporkov  et.  al, 
1998], 


Pmm=±- 


Cxx  (xm) 


+  <*»,)] 


Ax. 


The  upper  sign  corresponds  to  the  TM  case  and  the  lower  sign  to  the  TE  case  and 
^xx(xm)  is  surface  curvature  at  the  point. 


Direct  matrix  inversion  becomes  prohibitively  large,  requiring  the  storage  of  the 
NxN  propagator  matrix,  where  N  is  the  number  of  unknowns.  Furthermore,  the 
computation  time  for  LU  decomposition  scales  as  N3/3+N2  -5N/6  [Kapp  and  Brown 
1996];  here  decomposing  the  original  propagator  matrix  results  in  a  lower  triangular 
matrix,  L  and  an  upper  triangular  matrix,  U.  The  MOMI  approach  to  the  scattering 
problem  recasts  the  integral  equation  into  a  discretized  form  that  is  amenable  to  solution 
via  simple  forward  elimination  and  back  substitution  without  the  enormous  memory 
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requirements  of  LU  decomposition.  After  some  manipulation,  the  discretized  MFIE  can 
be  written  in  the  following  form 

J  =  [I  -  U  Y  [I  -  L]-1  V  +  [I  -  U]~!  [I  -  L]"1  LUJ  (8) 

Although  it  appears  that  matrix  inversion  is  still  needed  to  solve  (8),  it  can  be  shown  that 
alternating  forward  and  back  substitution  may  solve  this  equation.  The  first  term,  Jb,  has 
been  called  the  “new  Bom  term”.  The  following  is  a  general  iterative  solution  whose  first 
term  is  Jb  and  the  remaining  terms  are 

j = s  { &  -  urji  -  Lr  lu  }-[i  -  ur-fi  -  Lr1  j-  w 

n=0 

The  “new  Bom  term”  (n  =  0)  contains  all  orders  of  multiple  scattering  which  are 
continuously  forward  scattered,  continuously  backward  scattered,  and  those  which  are 
first  forward  scattered  and  then  backward  scattered.  Numerical  simulations  have  shown 
that  the  “new  Bom  term”  itself  is  adequate  for  most  practical  surfaces.  For  very  rough 
perfectly  conducting  surfaces,  a  maximum  of  two  MOMI  iterations  has  typically  proven 
to  be  sufficient. 
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2. 1  The  Incoherent  Power  via  the  Impulse  Response  Method 


The  use  of  the  impulse  response  method  has  been  well  established  in  literature  for  the 
calculation  of  the  average  incoherent  power  returned  from  the  ocean  surface  under  pulse 
illumination.  One  of  the  advantages  of  this  approach  is  numerical:  the  average  return 
power  can  be  recast  into  a  series  of  convolutions.  The  result,  consequently,  is  found 
easily  and  efficiently  using  the  Fast  Fourier  Transform,  the  FFT. 

The  power  due  to  an  incremental  area  with  a  given  backscattering  cross  section  is 
derived  directly  from  the  radar  equation.  Subsequently,  extending  this  power  to  include 
the  effects  of  the  entire  surface  will  lead  to  an  expression  for  the  returned  power  [Brown, 
1977].  The  average  power  returned  from  an  element  of  area,  dA,  with  a  cross-section  per 
unit  area  cr°(0,  <(>)  is  given  by  the  standard  radar  equation;  see  Figure  6. 


dP*(t)  = 


PT(t)7.2G2(0,(j>)g°(0,(t>) 

(4ti)3R4 


dA 


(1) 


where  R  =  range  from  the  radar  to  the  elemental  area,  dA 


G(0,  <j>)  =  antenna  gain  at  the  given  angles 
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The  average  power  returned  from  a  distributed  target,  such  as  the  terrain  in  this  case,  is 
calculated  by  a  superposition  of  backscattered  power  from  each  elemental  surface  area, 
dA.  A  superposition  of  power  is  appropriate  since  the  scattering  surface  is  assumed  to 
have  a  sufficiently  random  nature;  there  is  no  coherent  return.  The  average  backscattered 
power  returned  from  the  illuminated  surface  can  be  written  as 


i  2  »  2n 


Pft  2(ro  ~4(x>y) sec  8) 


P*(,)  ■  M  I!  V^ey  ■7°W(e.*)pd*dp  (2) 

where  the  slant  range  to  the  terrain  has  been  re-written:  R  =  r0  -  £(x,y)sec0.  Since  the 
integration  is  over  the  flat  surface,  cylindrical  coordinates  can  be  replaced  with  the  radar 
coordinates  (r02  =  p2  +  h2  =>  r0dr0  =  pdp),  and  the  range  to  the  surface  can  be 

approximated  by  the  mean  range  for  the  evaluation  of  amplitude  terms.  Hence,  the 
returned  power  can  be  re-written 

pr  L  _  2(ro  -£(x,y)sec  e)^ 

pR(t)  sA-JJ-i - ji! - iGI(e^)O0(e.*)r„#dr„  (3) 

W71/  0  0  ro 

The  average  power  returned  from  the  surface  can  then  be  found  by  evaluating  the 
expectation  with  respect  to  the  surface  heights. 
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2(ro-4(x,y)sec0) 


G!(0,'l>)7o(e>)ro#dropE©d5 


Introducing  the  change  of  variables, 


5(x,y)  = 


2  ^(x,  y)  sec  0 


P^  2sec0P^  2sec0^ 


we  express  the  average  returned  power  as 


i  2  ®  2 


f  r 


(P'(t)>  S  (4»J>  Hi 


PE(4)d?G2(0,<f>)o°(8,<t>)  r0d<j>dr0 


Recognizing  the  convolutional  form  in  the  random  variable,  \  ,  we  re-write  the  power 
received  as 


(Pr«) 


,  Pt 

*2  «  2  it  1 


M°l 

<N  | 

1 

®p? 

t— b 

1 

1 

l  co ; 

o 

O 

3 

*o 


G2  (0,  <j))a°  (6,  <j))d<|>  dr0 


(4) 


Substituting  for  the  constant  delay  term, 

t' 

and  introducing  the  shifting  properties  of  the  Dirac  Delta  Function,  the  power  received  is 
written  as  follows 


2ro 

Co 


A,2  %2*fPT(t-t’)®P?(t-t’) 


(PR(t))  =  J  - -G2(0, <j))a°(0, (j))d(j>dr0  5 

-ooV51/  0  0  r0 


f  2r  ^ 

f  _Z_£. 

V  co  J 


dt' 


”  I2  "r2?PT(t-t’  )®PF(t-t')  ,  , 

“  fzrvlF^ - - - G  (0.  (&.  4>)d<|> dr0  s 

iW  0  0  ro 


r  2r  ^ 

t'  — 

V  c0  J 


dt' 


Consequently,  the  average  scattered  intensity  from  a  rough  surface  can  be  expressed  as  a 
series  of  convolutions 

<Pr(t)>  =  PT(t)®p?(t)®PPS(t)  (5) 


From  equation  in  (5),  the  last  term  represents  the  average  backscattered  power  from  a 
transmitted  impulse  function  and  has  been  called  the  “Flat  Surface  Impulse  Response” 
(FSER) 


Pps(t)  = 


(4ti  y 


JJ 

Surface 


8|t-  — 


\ 

c  J 


G2(0,  tj>)a°(0,  <(>)dA 


(6) 
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where:  8(*)  = 

X 

G(9,  4)  = 

o(0,4)  = 

dA 

ro 

h 


a  Dirac  delta  function  which  accounts  for  the  two  way 
propagation  delay 

wavelength  of  the  carrier 

radar  antenna  gain 

surface  scattering  cross  section  per  unit  area 
elemental  surface  area,  dA  =  r0dr0d<j>  =  pdpd<j> 

slant  range  from  the  radar  to  the  mean  surface  at  dA 
radar  height  above  surface 
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3  Estimating  the  Return  Power  Components 

3.1  Introduction  to  the  use  of  a  Radiative  Transfer  Approach 


A  full  wave  approach  to  the  multiple  scattering  problem  presents  many  analytical 
and  numerical  challenges  even  in  tenuous,  i.e.  sparse,  media;  consequently,  the  simpler 
ideas  and  the  more  tractable  numerics  of  radiative  transfer  present  an  appealing 
alternative.  However,  since  phase  information  is  lost,  the  multiple  scattering  phenomena 
described  by  transfer  theory  are  not  well  understood.  In  addition,  transfer  theory  may 
only  have  a  certain  range  of  validity.  There  exist  at  least  two  different  levels  of  modeling 
the  environment  in  radiative  transfer  theory.  The  first  and  most  abundant  in  the  literature 
is  a  level  that  can  become  quite  detailed.  Typical  examples  of  this  are  found  in  references 
[Ulaby,  1990]  and  [Karam,  1997].  The  second  approach  treats  radiative  transfer  as  a 
theory  that  deals  in  bulk  media  and  effective  parameters.  A  typical  example  of  this 
approach  would  include  Schwering  [1985].  In  this  section,  we  present  the  formulation  of 
a  simple  radiative  transfer  model  for  predicting  the  return  power  waveform  from  a  rough 
surface  with  a  vegetated  cover  that  is  based  on  measured  parameters. 

A  computationally  efficient  method  for  the  determination  of  the  scattered  power 
density  can  be  found  in  [Brown,  1977]  for  a  rough  surface  and  [Newkirk  and  Brown, 
1996]  for  a  rough  surface  covering  a  penetrable  volume.  This  approach  creates  a 
numerically  efficient  method  since  the  incoherent  return  can  be  cast  as  a  series  of 
convolutions.  Henceforth  referred  to  as  the  Impulse  Response  approach,  this  method  as 
applied  to  rough  surfaces  has  been  briefly  reviewed  in  Chapter  2.  Although  the  Impulse 
Response  approach  was  originally  derived  from  the  radar  equation,  it  leads  to  a  volume 
return  waveform  that  can  be  derived  from  radiative  transfer  theory  [Adams  and  Brown, 
1998a].  This  idea  will  be  extended  to  the  volume  response  of  vegetation  over  a  rough 
surface  in  the  next  section. 
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3.2  Derivation  of  a  Volume/Surface  Impulse  Response  Approach 

The  geometry  for  the  radiative  transfer  approach  is  given  in  Figure  7  below.  In  this 
figure,  the  first  two  boundaries’  (enclosing  the  canopy)  height  statistics  are  described  by 
the  same  random  variable,  £(x)  and  the  third  (rough  terrain)  boundary’s  height  statistics 
are  described  by  the  random  variable,  £(x) .  The  mean  heights  of  the  layer  boundaries,  d] 
and  d2,  are  deterministic  distances.  Hence,  we  have  implicitly  assumed  a  zero  mean 
surface  with  a  layer  of  vegetation  whose  average  thickness  is  dj.  This  vegetative  layer  has 
mean  height  above  ground  equal  to  a  constant,  d2- 


Figure  7:  the  geometry  describing  the  volume  and  surface.  Note  that  dotted  lines 
indicate  average  levels  for  the  associated  boundary. 

Beginning  with  a  general  form  of  the  radiative  transfer  equation  for  the  incoherent 
power  density  or  the  intensity  in  the  medium,  a  simple  form  of  the  radiative  transfer 
equation  amenable  to  solution  via  convolution  will  be  derived.  The  geometry  for  the 
general  radiative  transfer  equation  is  given  in  Figure  8 
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Figure  8:  Scattering  geometry  for  the  intensity  [Ishimaru,  1997] 

Assuming  that  the  scattering  process  is  polarization  insensitive,  we  will  use  the  scalar 
radiative  transfer  equation,  which  relates  the  differential  change  in  the  power  density  over 
volume  ds.  This  is  written  as  (including  the  time  dependent  variation) 

=  -ptftI(s;f,t)  +  jj  p(s,s')l(s' ;  r,  t)dco' 

os  4n  4jc  ^ 

+  Js(s;f) - ^-I(s;f,t) 

cs(r)5t 

where 

•  I(  s ;  r )  is  the  power  density  in  the  s  direction  at  the  position:  r 

•  s  is  a  direction  of  the  power  density 

•  p  is  the  scatterer  density 

•  C7t(f)  is  the  scatterers  total  cross  section  which  is  the  sum  of  the  absorbing  and 
scattering  cross  sections:  at(r)=  crabs(r)  +  crsc(f)and  as  written  here,  may  be  a 
function  of  position  r . 

•  p(s,s')  is  the  scattering  function  of  each  scatterer;  (prime  denotes  incident 
direction(s))  and  is  related  to  the  amplitude  of  the  field  scattering  function  squared. 

•  Js(s;r)  is  the  source  function  (emission  sources) 
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Referring  to  equation  (1),  the  change  in  power  in  the  r  direction  is  proportional  to  the 
power  incident  on  the  differential  volume  element.  This  power  is  then  depleted  by 
absorption  as  well  as  scattering  into  other  directions.  On  the  other  hand,  the  power,  as  it 
propagates  through  the  differential  volume,  increases  by  an  amount  due  to  scattering  into 
the  r  direction  from  other  directions  r'  as  well  as  energy  emitted  inside  the  differential 
volume: 

We  can  now  derive  an  impulse  response  representation  by  making  the  following 
assumption  regarding  the  scattering  function  (or  classically,  the  phase  function),  p(r ,  r ’)  . 
It  will  be  assumed  that  each  scatterer  scatters  energy  in  the  forward  and  backward 
directions,  exclusively. 

p(s,s')  =  —  [af5(s'-s)  +  cb8(s'  +  s)] 
a. 


where  af  and  Gb  are  the  position  dependent  forward  and  backward  scattering  cross  section 
of  each  scatterer,  respectively.  These  may  be  functions  of  depth  into  the  media,  shown 
explicitly  by  the  r  dependence,  as  well  as  the  scattering  angle.  In  addition,  we  will 
assume  that  there  are  no  emission  sources  present;  consequently,  the  source  term, 
Js(s;r),  is  zero.  This  is  a  good  assumption  for  active  sensing  techniques  [Ulaby,  1986]. 

With  this  scattering  function  assumed,  the  radiative  transfer  equation  is  recast  into  a 
greatly  simplified  form.  Since  the  direction  of  power  density  propagation  s  has  been 
limited  to  the  radial  direction,  r ,  the  equation  governing  the  power  density  can  be  written 
as  a  first  order  partial  differential  equation  in  two  variables:  time  and  distance.  Implicitly 
assuming  the  r  dependence  in  the  cross  section  parameters,  the  simplified  equation  of 
transfer  becomes 


ai(r;r,t) 

dr 


-pcrtI(r;r,t)  +  p[afI(r;f,t)  +  abI(-r;r,t)]  - 


1  dl(r,r,t) 

cs(f)  at 


In  order  to  further  simplify  the  equation  of  transfer,  we  split  into  it  into  two  parts: 
downwelling,  that  power  density  which  propagates  in  the  forward  hemisphere  and 
upwelling,  that  power  density  which  propagates  in  the  backward  hemisphere  as  defined 
by  the  direction  of  propagation,  r . 
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Let  us  first  consider  the  downwelling  intensity.  In  its  solution,  we  will  assume  that 
the  upwelling  power  density  does  not  act  as  a  source  for  the  downwelling  power  density 
or  ab(r)  =  0 .  At  this  time,  we  consider  this  event  to  be  second  order  scattering  that  may 
be  neglected.  Next,  we  define  an  effective  extinction  coefficient  per  unit  volume 
ke(r)  =  pcrt(r)  -paf(r) .  We  then  assume  a  media  that  is  radially  distributed  which 

results  in  the  modified  effective  extinction  coefficient:  ke(r)  =  p<rt(r)  -paf(r).  Hence, 
implementing  these  assumptions,  we  find  the  greatly  simplified  equation 


dl(f;r,t) 

dr 


(r)I(f';r,  t) 


1  dl(f;f,t) 
c.(?)  at 


(2) 


Since  it  has  been  assumed  that  the  upwelling  power  density  does  not  contribute  to  the 
downwelling,  there  is  no  coupling  of  power  from  the  upwelling  into  the  downwelling. 
Consequently,  given  an  initial  power  density  at  the  upper  foliage  boundary  with  free 
space,  I0(t  -  t’)  where  t’  is  a  range  dependent  delay,  the  solution  for  the  downwelling 
power  density  is  found  in  closed  form.  The  method  of  characteristics  yields  a  time-shifted 
argument  for  the  power  density,  while  the  distance  dependence  can  be  found  by  simple 
integration. 


I(r ;  r,  0,  <J>,  t)  =  Ic 


t- 


(3) 


where  Io  (t  - 1’)  is  the  time-delayed  incident  power  envelope;  the  time  delay  is  a  function 
of  the  range  in  free  space  from  the  antenna  to  the  canopy  (ri)  and  the  range  into  the 
medium  which  may  have  a  range  dependent  group  velocity,  cs(r).  Note  that  the 
downwelling  power  density  is  directed  in  the  r  direction  or  along  a  radial  path  from  the 
source  antenna. 


The  differential  equation  governing  the  upwelling  power  density  has  a  similar  form; 
however,  the  downwelling  power  density  acts  as  a  source  for  the  upwelling.  In  addition, 
the  upwelling  power  density  is  directed  in  the  -  r  direction  or  along  a  reverse  radial  path 
toward  the  source  antenna.  Consequently,  the  governing  differential  equation  is  the  same 
with  the  exception  of  the  coupling  term  relating  the  upwelling  and  the  downwelling 
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intensities.  The  differential  equation  governing  the  upwelling  power  density  is  expressed 
below. 


dl(-f;r,t)  _  r  ,  *T,  f  . _ 1  dl(-r,s,t)  _  f 

Q  Ke(r)I(  r,r,t)  +Gb(r, 0, (p)I(r, r,t) 

dr  cs(r)  at 


Subsequently,  substituting  the  solution  for  the  downwelling  power  density  into  equation 
(4),  the  following  differential  equation  is  created  governing  the  upwelling  power  density 


-  -k,(r)I(-r;r,«)  - 

ar  cs(r)  at 


+  CTb(r,e,<j>)I0  t-  — -  r-^-)exp{-  f  ke(p)dp} 

v  co  ,csWJ  ' 


The  attenuated  downwelling  power  density  that  passes  through  the  foliage  layer  and  is 
subsequently  scattered  by  the  underlying  surface  acts  as  a  source  for  the  upwelling  power 
density  at  the  foliage  layer’s  lower  boundary.  In  addition,  the  downwelling  power  density 
continuously  contributes  to  the  upwelling  power  density  due  to  the  coupling  term.  Note 
that  this  source  was  absent  in  the  differential  equation  for  the  downwelling. 
Consequently,  in  formulating  the  solution,  the  upwelling  power  density  has  two 
independent  sources:  the  power  waveform  backscattered  by  the  surface  and  the 
backscattered  downwelling  power  density  from  within  the  volume.  Finally,  the  upwelling 
power  density  is  evaluated  at  the  top  of  the  canopy  (r  =  ri).  Again,  invoking  the  method 
of  characteristics  as  a  solution  method  for  the  time  dependence,  the  resulting  solution  of 
equation  (6)  has  two  independent  terms 


where 


as(0,<j>)  =  a°(0,<j>)dA,  the  surface’s  scattering  cross  section 
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ri  =Iio+4(x)sec0 

r2  =  r10  +  d,  sec  0  +  £(x)  sec  0 

r3  =  rI0  +  (d,  +  d2  )sec  0  +  £(x)  sec  0 

This  expression  shows  a  simple  superposition  of  two  terms;  the  first  term  is  the  rough 
surface  return  propagated  back  up  through  the  foliage  and  the  second  term  represents  the 
foliage  scattered  return.  In  order  to  construct  the  average  return  power  in  the  impulse 
response  format,  these  two  responses  are  averaged  and  manipulated  to  yield  an  impulse 
response  term  in  each  case.  However,  an  additional  assumption  is  necessary  for  a  fully 
convolutional  result  similar  to  that  given  in  the  literature  for  a  rough  surface  alone:  the 
random  variables,  <^(x)  and  c(x) ,  describing  the  canopy  and  the  rough  surface, 
respectively,  must  be  assumed  to  be  independent. 
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3.3  Incoherent  Scattered  Power:  the  Volume  (Foliage)  Return 

In  the  formulation  of  the  scattering  from  a  rough  surface  with  a  vegetative  cover,  we 
have  assumed  that  scattering  occurs  exclusively  in  the  forward  and  backward  directions; 
this  implied  that  the  power  density  in  radial  direction  r  does  not  interact  with  the  power 
density  in  any  other  radial  direction.  This  in  turn  has  led  to  a  closed  form  result  for  the 
downwelling  intensity  and  consequently,  an  uncoupled  relatively  simple  equation  for  the 
upwelling  power  density,  equation  (6)  in  the  previous  section.  The  two  terms  of  the 
solution  in  (6)  can  be  simplified  independently.  Each  represents  a  different  scattering 
phenomena,  surface  and  volume  scattering.  In  this  section  we  examine  the  foliage  or 
volume  return. 

We  begin  with  the  second  term  of  equation  (6)  for  the  upwelling  power  density,  the 
volume  response.  After  substitution  for  the  slant  range  variables  (rj,  T2,  ...)  with  the 
associated  distance  and  random  variables  as  a  function  of  antenna  pointing  angle,  0, 

ri  =rio  +  4(x)sec0 

r2  =  r10  +  d,  sec  0  +  §(x)  sec  0 

r3=r10  +  (d1  +d2)sec0+  £(x)sec0 


the  power  density  is  found  to  be  approximated  by  the  following 


fT,o+E(x)sec0+d,sec6  )  i-a 

I(~f;r,0,<t>,t)  =  f  B  ab(a)exp  -2f 

Jrm  +  5(x)sec0  '  Jr io 


■|O+5(!Osec0i 


ke(n)dp 


.  J  (t  (rio  +  ^(x)sec  0)  _  2  p 


5(x)sec0  c(p) 


In  general,  following  a  slightly  modified  version  of  the  method  of  Adams  and  Brown 
[1998a],  and  assuming  a  layered  media  with  parallel  boundaries,  the  average  power 
density  can  be  expressed  as 


(  (ri0  +^(x)sec  0) 
°  <  c0 


dap^(4)d^ 


(2) 
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In  order  to  create  a  convolutional  form,  the  integration  limits  must  be  extended  to 
infinity.  The  upper  limit  may  be  extended  to  infinity  by  assuming  that  the  extinction 
coefficient  becomes  very  large  once  the  range  extends  beyond  the  lower  foliage 
boundary,  this  will  effectively  eliminate  the  volume  return  after  the  lower  foliage 
boundary  is  surpassed.  The  lower  limit  of  integration,  on  the  other  hand,  can  be  extended 
by  the  use  of  the  unit  step  function,  u(r  -  rj).  Consequently,  the  average  power  density  can 
be  rewritten  in  terms  of  integrals  with  infinite  limits.  Under  the  change  of  variables, 
p'=  p  -  [r10  +  ^(x)  sec  0],  the  expression  for  the  upwelling  power  density  becomes 


(l(— ?;r,  9, <(>, t))  =  J  Jab(a)I0 


t  (rio  +  4(x)sec6)  _  2  ra-[r10+S(x)sece]  dp’ 
c0  0  c(p') 


•exp{-2j0a  Iri°+^(x)sec0]  ke(p')dp’  }u(a-  [r10  +£(x)sec0])dap  ^(£)d^ 

(3) 

Assuming  that  there  is  no  volume  return  from  the  atmosphere  between  the  antenna  and 
the  foliage  crown,  the  backscattering  cross  section,  which  is  a  function  of  distance,  can 
also  be  shifted  by  the  slant  range.  Defining  two  new  functions 


(4) 


E(Y)  =  ab(y)exp|  -2  Jke(p')dp’ |  u(y ) 


(5) 


the  average  upwelling  power  density  at  the  upper  foliage  layer  can  be  rewritten 

f  \ 

t  - 

-00  -00  V  Cr 


(l(— f;  r,0,  <j>,t))  =  Jlljt-i-  g[a  -  ri  ]  E(a  -  r,  )da  p  %  (£)  d^ 


(6) 


where  it  has  been  previously  defined  that  r,  =  (r10  +^(x)sec  0).  Following  the  method 


of  Adams  and  Brown,  [1998a],  the  following  definitions  are  constructed  which  transform 
distance  into  time 
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t'  =  gja  -  r,  ] 


=>  a-r,  =  g-'(t') 


dt’=  g’[a-r,]da  =  g,(g‘1(t’))da 


Substituting  these  expressions  into  the  average  upwelling  power  density  of  equation  (6), 
the  average  upwelling  power  density  is  reconstructed  in  the  following  form 

<I(-f;r,e,f  t)>  =  j  j  lit  -i  -  t'|^ff  ’^dt'P;(0<lS 

-CO  -00  l  C0  /  g  ^g  )) 


=  J  J  I0  t  — i E(t')dt'p,(5)d^ 


where  a  new  function  has  been  defined: 


Bo  =  MM 

g'(g-'(t)) 

Noting  that  expression  (7)  contains  a  convolution  in  the  variable  z,  we  perform  the  z 
integration,  leaving  the  result  in  the  form  of  a  convolution  (with  convolution  represented 
by  the  symbol:  <8>)  shown  in  brackets  below 

.  f  (  r  ^  J  r  Yl 

(I(-f;r=r„e,<M))=  /  I0  t ®E  t  -- Us©d| 

’*  i  V  c0  J  V  coJ J 

Substituting  for  the  slant  range  in  terms  of  the  distance  to  the  mean  height  and  the 
random  variable  representing  the  distribution  about  the  mean,  i.e.  r,  =  (r,0  +  £(x)sec0) 


<I(-f;r  =  r„e,*,.»=  ?  jl/ .  W . 

-[  l  C0  J  { 


+  £(x)sec0) 


First,  we  substitute  for  the  constant  delay  term:  we  let 
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Then  we  make  a  change  of  variables  with  respect  to  the  random  variable  representing  the 
crown  height  statistics;  we  form  a  new  random  variable  and  its  associated  probability 
density  function 

^  £(x) sec  0 

5(x)  =  — -  , 


P|(t)  = 


"_£oLl 

^secGj 


Consequently,  the  average  upwelling  power  density  becomes 

(I(-f;r=r„e,+,t))  =  -&-]  {l„(«-t,-5  )®E(t -t„- ?)}p,(I)d| 

SCC  t7  -oo 


Again,  the  average  upwelling  power  density  is  re-expressed  in  the  following 
convolutional  form  with  respect  to  the  modified  surface  height  random  variable 


(l(-f;r  =  r„0,*,t))  =  -\l0(t -t„  )®E(t -t0)  (8) 

sec  0  ? 

Finally,  in  order  to  find  the  total  power  returning  toward  the  radar,  we  integrate  the  power 
density  over  a  surface.  In  this  case,  a  convenient  surface  is  the  top  of  the  canopy. 
Allowing  for  full  penetration  of  the  incident  power  density  (i.e.  no  reflection  from  the 
boundary  separating  free  space  from  the  foliage),  we  substitute  for  the  incident  power 
density  at  the  canopy  (expressed  in  the  {  }  brackets  below)  weighted  by  the  antenna  gain 
G(0,  <)>)  in  the  direction  (0,  (()) , 

47ir10 

Furthermore,  we  must  assume  a  narrow  beamwidth  such  that  sec  0  »  sec  0O  (boresight 
direction:  0o).  This  is  required  since  the  integral  to  be  performed  over  the  radial 
coordinate  implicitly  contains  0  dependence;  otherwise,  a  convolutional  form  can  not  be 
obtained.  Allowing  for  the  additional  delay  due  to  the  transmission  back  to  the  antenna 
(an  additional  to)  from  the  canopy  and  the  receiving  antenna’s  effective  aperture,  the 
following  result  is  obtained 
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(p«<‘»  =  ]/^V-j^rrlpT(t-2to)}®E(,-2t0)®p5(t-2t0) 


0  0  sec9o  l  4™io 


X2G(0,<j>)  jaj 
■/-  "fr  rio  d^10 


substituting  the  original  expression  for  to,  we  find 


(PR(t)>  = 


0(0,4.) 


t--=32-  ®E  t-^  ®pJt-^fi.  r10 d<jxir.0 

secG0  o  o  1  47crI0  (  c0  Jj  (  c0  J  c0  J  (4rt)  r2 


2r10  U2G(0,<{>) 


=  f/v  ‘a  H j^riMlp  f t 

(47t)  sec  0O  o  o  [  r10  (  c,  ))  (  c0  )  « 


2r  i 

t-±iLi  d4>dr1( 


Exploiting  the  shifting  properties  of  the  Dirac  delta  function  and  rearranging  the  resulting 
integrals  yields 

jy^6r^Wdt. 

o  o  r.n  l  1 


=  J jpT(t  -t')®E(t  -t')®ps(t  -f)}pfs(t')  dt' 


-  PT(t  )®E(t  )®p^(t  )®PFs.(t) 

where  the  transmitted  power  waveform  is  given  by  Px(t)  and  the  modified  Flat  Surface 
Impulse  Response  function  (see  Chapter  2)  is  given  by 

p's'(t)=^^*’S^4,~~]d*dr“  (9) 

(47i)  sec  0  0  0  r10  (  c 0  ) 

the  modified  probability  density  function  for  the  crown  height  statistics  is  given  by 

P*(0  =  0°) 

*  sec0o  (sec  B0  J 
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and  E(t)  is  a  function  relating  decay  to  depth  of  penetration  into  the  foliage  layer 


E(t)  = 


Efe^t)) 

g'(g_1(t)) 


(11) 


This  general  solution  for  the  volume  response  can  be  modified  to  yield  a  simpler 
result.  Assuming  that  the  velocity  is  constant  in  each  layer  of  the  medium,  the  solution 
becomes  more  apparent.  Here,  we  assume  that  medium  1  contains  the  leaves  and 
branches  (group  velocity  is  cvi)  and  medium  2  is  the  trunk  region  (group  velocity  cV2)- 
The  starting  point  for  the  upwelling  power  density  due  to  the  volume  return  (along  a 
radial  in  the  r  direction)  is  then  given  by 


_ .  A  „  .  .  frlo+5(x)sec0+d,sece 

I(— r;r,e,«t>,t)  =  [  ^(a)1 

Jr,„+E(x)sec8, 


t  - 


(r10  +^(x)sec6)  l[a  -  (r10  +  %(x)  sec  9)] 


'vl 


•expj-2|  k.(u)du  (da 

r  t  Jr,0+4(x)sece,  e  n  ) 

After  following  the  previous  procedure,  the  average  power  as  a  function  of  time, 
scattered  from  a  volume  with  an  irregular  interface  at  the  crown  can  be  expressed  in  the 
convolutional  form 

(PM)  =  PT  (t)  ®  Pre(t)  ®  P?(t)  ®  Eft) 


where  in  this  particular  case, 

t'  =  g(r')  =  —  =>  r'  =  g_1(t')  =  —  and  g'(r)  =  —  =  constant 
c  2  c 


Hence,  for  a  group  velocity  in  the  volume  given  by  cv, 


E(t)  = 


E(g-’(t))  E(cvt/2) 

g'(g“l(t))  2/cv 


y  ob  exp{  “  2!o't,Z  ke(p') dp' }  u(cvt  / 2 ) 


Note  that  the  unit  step  function  u(t)  “turns  on”  when  t  =  0. 
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3.1  Incoherent  Scattered  Power:  the  Rough  Surface  Return 

We  start  with  the  expression  for  the  power  density  attenuated  by  the  foliage  in  its 
downward  passage,  scattered  from  the  surface,  and  then  attenuated  by  the  foliage  in  its 
upward  passage;  this  is  the  first  term  of  equation  (6)  in  Section  3.1.  Note  that  the 
geometry  of  Figure  3.1-1  still  applies 

I(-?;r  =  r„e,<M)  =Cfs(G, <(►)!<,  t_~_2f~7T  exp{“2  Pke(p)dp  }  (1) 

<  c0  1  C(W  y  1 

Substituting  for  the  power  density  using  the  following  relationship 


I(r,0,<M) 


PT(t)  G(9,4>) 
4nr2 


(2) 


where  I  is  the  power  density  (sometimes  called  intensity),  PT  is  the  power  waveform  and 
the  r  direction  is  specified  by  the  angles  0  and  <j>.  The  average  power  returned  from 
within  the  illuminated  region  can  be  evaluated  by  integrating  over  a  surface 
encompassing  the  illuminated  area.  In  this  case,  we  choose  to  integrate  over  the  area  at 
the  top  of  the  canopy  (r  =  ri).  Hence,  substituting  the  power  waveform  for  the  incident 
power  density  in  equation  (1)  via  the  relationship  in  equation  (2)  and  performing  the 
ensemble  average  over  the  random  variables,  the  total  average  power  at  the  crown  is 


(3) 

where  the  surface  scattering  cross  section  per  unit  area,  (0, 4>) ,  has  been  included.  The 

angles  (0,4>)  are  spherical  coordinates  centered  at  the  antenna  and  can  be  related  to  the 
variables  of  integration.  Also  in  this  expression,  the  antenna  gain  has  a  boresight  angle 
given  by  (0O,4»O)  and  the  angles  (0a,4>a)  are  spherical  coordinates  defined  with  respect 
to  this  antenna  boresight  direction  which  may  also  be  related  to  the  variables  of 
integration. 
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Although  the  solution  procedure  can  proceed  with  a  propagation  speed,  c(r) ,  and  an 
effective  extinction  coefficient,  ke(r)  that  vary  with  position  as  assumed  in  (3),  the 
following  results  are  simplified  since  they  are  based  on  a  constant  group  velocity  and 
extinction  coefficient  in  each  layer. 

•  Free  space,  the  group  velocity  is  co 

•  Layer  1,  the  canopy  region,  the  group  velocity  is  cvi,  the  effective  extinction 
coefficient  is  kel,  r  6  (r,,r2) 

•  Layer  2  the  trunk  region,  the  group  velocity  is  Cv2,  the  effective  extinction  coefficient 
is  ke2,  r  e  (r2,r3) 


Hence,  the  integrals  with  respect  to  the  radial  distance  within  the  argument  of  the 
transmitted  power  may  be  easily  performed,  yielding  the  average  power  at  the  radar 


(p(t»  = 


jj  K<e,« 


surface  at 
r  =  rio 


A,2G2(ea,(j)a) 

(47t)3r4 


PT 


2ri  2(r2  ~  ri )  Zfc-rz) 


■'vl 


'v2 


exp{- 2  r 


ke(p)dp  dS 


The  final  results  will  be  cast  in  a  convolutional  form  for  the  average  returned  power. 
After  expanding  the  transmitted  power  waveform’s  delay  time  argument  in  terms  of  the 
random  variables  and  constant  terms, 

ri  =rio  +^(x)sec® 

r2  =r10  +  d,  sec0  +  i;(x)sec0 

r3  =  r10  +  (d2  +d2)sec0+  ^(x)sec0 

and  performing  the  integrations  with  respect  to  the  extinction  coefficients,  the  average 
return  power  as  a  function  of  time  becomes 
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/n  /*\\  U  >-:G:f6;,<pa)  c 

JJ  7233T  g.(9''tl> 

Surface  at  V*KJ  1 

f  =  f|0 

r°  p  p  t  2rlosec0  £(x)sec0  2d,  sec  0  2d2sec0  2[^(x)-£(x)]sec0 

J-co  J-oc  ^  p  p  r  r  p 


exp{-2keld,  sec0-2ke2(d2  +  £(x)-£(x))sec0  ]p5(^)p?(Qd£d<^rdrd<j> 


where  p,(4)  andpc(Q  are  the  probability  density  functions  for  the  boundary  heights  of 

the  foliage  volume  and  the  rough  surface,  respectively.  Rearranging  and  substituting  for 
the  constant  terms,  equation  (4)  is  rewritten 

<PR  ('» =  JJ  ^J%#a»(e,1|1)exp{-2sece(kI1d1-k.,d2)} 

Surface  at  V  ^ 

r=rl0 


££p. 


2£,(x)sec0  2i^(x)sec0 


■exp{2kc24(x)sec0  }p^)d£exp{-2ke2C(x)sec0  }pc(Qd;rdrd<|> 
where  ke,,ke2are  the  effective  extinction  coefficients  in  medium  1  and  2,  respectively; 


_  2r,osec0  2d,  sec  0  2d2sec@ 


2  (  1  2 


Ca  VC0  Cv2 


Performing  a  change  of  variables  which  transform  distance  into  time 


2^(x)sec0 


2^(x)sec0 


—  v2 


2sec0 


d£,  =  i —  dt2 

2sec0 


and  defining  some  new  functions, 
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ftiOi)  =  exp{-ke2cv2t1}ptl(t1) 
ft2(t2)  =  exp{  ke2  cat2  }pt2(t2) 


Noting  that  the  probability  density  functions  describing  the  boundaries  must  also  be 
transformed, 


Pti(ti)  = 


Pt2(t2)  = 


2sec0 


we  find  that  the  average  returned  power  can  be  expressed  as 

(Pr(0)=  JJ  ^  % T1 * '  Gs (0» 40 exp {-  2 sec 0(k£A  - ke2d2 ) } 

Surface  at  4W  ri0 

r=rio  (6) 

00  00 

T1  o  J  JPt(^  —  —  tj  — 12  )ftl(tI)dt1  ft2  (t2)dt2  r10  dr10d<j) 

4  sec  o  ^  ^ 

Performing  the  ti  and  the  t2  integrations  and  expressing  the  result  in  convolutional  form, 
which  is  represented  with  the  ®  symbol  we  find 


ft(*))=  JJ 


Surface  at 


cacv2  A,  G  (0a,(j)a) 
4  sec 2  0  (4Tt)3r140 


a°  (0,  ()>)  •  exp{-  2  sec  0(keld,  +  ke2d2 )} 


r  =  r 


10 


•  [PT  (t  - 10 )  ®  f«  (t  - 10 )  ®  ft2  (t  - 10)]  r10  dr10  d<j) 

Integrating  over  the  foliage  upper  surface,  the  returned  power  can  be  recast  into  the 
following  convolutional  form,  following  the  methods  outlined  in  previous  section. 

<P,(t)> = PT  (t)  ®  f„  (t)  ®  f,2  (t)  ®  PK.(t)  (7) 

where  ftl  (t)  and  ft2  (t)  are  functions  which  depend  on  the  probability  density  functions 
whose  random  variables  are  functions  of  the  random  variables  representing  the  surface 
and  canopy  statistics,  £(x),  <^(x)  as  well  as  the  extinction  coefficients,  and  the  antenna 
boresight  angle,  0.  The  flat  surface  impulse  response  function  (FSIR),  Pfs’O),  is  similar  to 
the  standard  FSIR  with  the  modifications  (among  others)  that  account  for  attenuation: 
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(8) 


p  m  yf8(’-2r%)G!(9„,»,)CT;(e,^) 

FS  '  4(4^  J  J  (r')4  sec1 0 

•exp{-2sec0(ktld1  +ke2d2)}r'dr’d(f) 

Here  the  antenna  gain  is  approximated  by  a  circularly  symmetric  pattern  with  a  pointing 
angle  given  by  (0o,<j)o)  and  the  angles  (0a,<t>a)  are  spherical  coordinates  defined  with 
respect  to  the  antenna  boresight  direction.  Consequently,  the  antenna  gain  can  be 
represented  by 

G(0a,<t>a)  =  Go(0o,<Mexpj-^sin20aj  (9) 

This  expression  (8)  for  the  FSIR  includes  additional  effective  parameters  related  to  the 
speeds  in  the  different  media,  the  incidence  angle,  and  consequently,  time  or  range: 

h'  =  — d,  +  — d2  +h 

Cvl  Cv2 


(r')2  =  (h')2  +  p2 

The  integral  for  the  FSIR  can  be  simplified  using  the  method  presented  in  Brown  [1977]; 
we  begin  by  substituting  the  two-way  incremental  ranging  time  for  the  actual  time: 

t=  t  ~^yc  ■  Assuming  that  the  beam  is  narrow  such  that  the  surface  incremental  cross 

section  is  constant  over  the  angular  extent  of  interest  and  that  sec0  =  sec0o  (boresight), 
the  FSIR  is  found  in  the  following  form 


?.2cocacv2(h')2Go(9o,<|>o)Cs(0o,<t>o) 

rFS’Vxl  ~  77  75  exP 

167t3(c0x  +  2h')5 


c0x  +  2h' 


(keidi  +ke2d2) 


In  4 

J  exp(  — sin2  0a  )d<() 

o  l  r  J 


Note  this  expression  can  include  an  asymmetrical  antenna  pattern  [Newkirk  and  Brown, 
1992], 
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Summarizing,  the  final  form  of  the  solution  to  the  average  power  returned  due  only 
to  the  rough  surface  can  be  expressed  in  the  following  convolutional  form, 

(p,  (t)> = pT  (t)  ®  f„  (t, )  ®  4  (t2 )  ®  pra.(t)  (ii) 

where  ftl(t,)and  ft2(t2)  are  functions  of  time  which  depend  on  random  variables  which 
are  functions  of  the  surface  and  canopy  statistics,  £(x),  ^(x) ,  as  well  as  the  extinction 
coefficients,  and  the  pointing  angle: 


for  layered  media  with  a  constant  velocity  in  each  media  and  where  ke],  ke2  are  the 

effective  extinction  coefficients  in  medium  1  (foliage  region)  and  medium  2  (trunk 
region),  respectively. 
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4  Interaction  between  the  Foliage  and  the  Surface 


The  returned  power  from  the  volume  must  be  modified  by  the  addition  of  the  time- 
delayed  return  from  the  rough  surface  under  the  random  media.  As  previously  seen,  the 
radiative  transfer  method  accounts  for  this  second  scattering  event  simply  by  assuming 
that  the  incident  wave  to  the  rough  surface  is  due  to  an  attenuated  version  of  the  original, 
free-space,  time-delayed  incident  power  waveform.  The  attenuation  is  due  to  the 
collection  of  scatterers  along  each  radial  from  the  antenna.  After  the  power  waveform  is 
then  scattered  by  the  surface,  it  again  travels  back  up  through  the  foliage,  suffering  the 
same  attenuation.  The  details  of  this  approach  were  developed  in  the  Chapter  3.  A  sample 
set  of  volume  and  surface  scattered  waveforms  based  on  the  radiative  transfer  results  of 
Chapter  3  is  given  in  the  first  section  of  this  chapter. 

Since  many  assumptions  are  inherent  in  the  radiative  transfer  result,  there  is  a 
question  as  to  the  validity  of  this  approach.  Can  the  interaction  between  the  foliage  and 
surface  be  modeled  simply  by  this  single  interaction?  The  foliage  scatters  a  field,  Ef , 
toward  the  antenna  and  other  directions.  This  foliage  scattered  field,  Esf ,  is  also  incident 

on  the  rough  surface  in  addition  to  the  free  space  incident  field  E,nc  (see  Figure  9).  It  is 
then  scattered  back  through  the  foliage  (see  Figure  10)  resulting  in  a  second  scattered 
field  returned  to  the  radar,  Esfsf ,  due  to  the  foliage-surface-foliage  interaction. 


Es 


ct.  ^ 


Figure  9:  The  total  incident  field  with  respect  to  the  surface 
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Naturally,  the  question  arises:  are  the  multiple  interactions  between  the  foliage  and 
the  surface  a  necessary  component  in  this  portion  of  the  model?  That  is,  is  there  a 

significant  second  order  surface  interaction  due  to  incidence  of  the  field  E^f ,  the  foliage- 
surface-foliage  field,  to  the  surface?  In  other  words,  when  E4f  is  incident  on  the  surface 
(note  the  addition  of  Efsf  into  Figure  11  with  respect  to  Figure  10)  will  there  be  a 
significant  correction  to  the  surface  scattered  field?  This  additional  incident  field  will 
modify  the  surface  currents,  which  will  radiate,  Efsfs ,  creating  a  new  foliage  scattered 

field,  EfSfsf .  The  third  order  approximation  to  the  interaction  between  the  foliage  and  the 

surface  would  repeat  this  process  again.  This  process  will  continue  indefinitely,  or  until 
the  corrections  become  negligible.  Notice  that  with  each  iteration,  the  final  foliage 
scattered  field  is  not  used  as  an  incident  field  for  the  surface. 
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Figure  1 1 :  the  second  order  approximate  scattered  field  from  the  foliage  and  surface 

combination 


Hence,  extensions  of  this  single  passage  event  would  include  an  infinite  series  of 
interactions  between  the  foliage  and  the  surface  in  the  case  of  single  scattering  theory  or 
fully  coupled  integral  equations  in  the  integral  equation  approach.  Is  the  first  order 
scattering  interaction  term  adequate?  Three  methods  of  approximation  to  the  second  term 
have  been  or  will  be  examined  in  this  study: 


1 .  Modified  “First  Order  Multiple  Scattering”  theory 

2.  Exact  formulation  using  the  Method  of  Moments  (MOM) 

3.  Reduced  integral  equation  approach 

Due  to  the  complexity  of  the  second  approach  listed  only  a  limited  number  of 
scatterers  can  be  placed  above  the  surface.  Hence,  we  will  investigate  the  higher  order 
interactions  based  on  a  single  scatterer  above  a  rough  surface.  Since  the  radiative  transfer 
approach  cannot  simulate  this  situation,  another  approach  was  required.  After  some 
investigation  the  first  listed  approach  was  found  to  not  only  support  the  single  scatterer 
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investigation,  but  also  turned  out  to  be  a  more  general  version  of  the  radiative  transfer 
result  as  developed  in  Chapter  3.  This  first  order  multiple  scattering  approach,  under 
certain  assumptions  reduces  to  the  result  given  for  the  radiative  transfer  approach.  This  is 
detailed  in  Appendix  A. 

The  first  approach,  the  modified  first  order  multiple  scattering  result,  does  begin  with 
a  single  scatterer.  Hence,  since  the  convolutional,  radiative  transfer  approach  is  related  to 
this  method,  we  need  only  show  that  the  foliage  to  surface  to  foliage  interaction  requires 
only  the  first  order  interaction,  i.e.  truncate  the  infinite  series  of  interactions  as  previously 
described  with  only  the  first  interaction  to  verify  the  assumption.  A  proposed  method  for 
extending  this  approach  to  a  collection  of  discrete  scatterers  is  also  outlined.  In  addition, 
if  this  method  were  successfully  implemented  as  a  convolution,  it  would  serve  as  a  more 
general  approach  than  the  radiative  transfer  method  as  developed  in  the  first  section  of 
this  chapter. 

Verification  of  the  first  order  multiple  scattering  result  will  require  a  comparison  with 
an  exact  solution.  Consequently,  the  next  section  of  this  chapter  examines  the  exact 
formulation  for  a  single  scatterer  above  a  rough  surface  and  solution  via  the  efficient 
MOMI  method  as  previously  described.  This  result  may  serve  as  an  exact  result  when 
compared  with  the  first  order  multiple  scattering  solution  obtained  for  a  single  scatterer 
above  a  rough  surface. 

Finally,  in  the  following  section,  the  exact  integral  equation  method  is  simplified 
using  some  reasonable  assumptions.  This  method  will  result  in  a  more  accurate  method  to 
simulate  interaction  between  a  single  scatterer  above  a  rough  surface  than  the  first  order 
multiple  scattering  approach.  In  addition,  it  may  also  yield  a  more  tractable  numerical 
model  when  it  is  extended  to  a  collection  of  scatterers  above  a  rough  surface  than  the  full 
Method  of  Moments  approach.  This  more  accurate  representation  of  the  interaction  may 
be  required  at  some  level  of  simulation. 
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4.1  Capabilities,  Limitations  and  an  Example  Return  Waveform 
for  the  Radiative  Transfer  Approach 

Previously  noted  limitations  of  the  radiative  transfer  result  have  included  a  limited 
scattering  pattern  for  each  volume  scatterer  and  a  narrow  beamwidth  approximation.  The 
chosen  scattering  pattern  demands  strictly  forward  scattering  and  backward  scattering. 
This  assumptions  decreases  the  number  of  coupled  differential  equations  from  N  (when  N 
scattering  directions  are  used  in  a  quadrature  approximation  to  the  integral  in  the 
radiative  transfer  equation)  down  to  two:  coupled  integral  equations,  one  governing 
forward  and  one  reverse  scattering.  Secondly,  the  upwelling  power  density  is  assumed 
not  to  influence  the  downwelling.  This  last  assumption  is  key  since  it  allows  a  closed 
form  solution  for  the  downwelling  power  density.  Otherwise,  the  solution  would  be  in  the 
form  of  two  coupled  differential  equations. 

In  Appendix  A,  the  radiative  transfer  result  is  shown  to  be  equivalent  to  the  first  order 
multiple  scattering  result.  From  the  first  order  multiple  scattering  analysis  of  Appendix 
A,  another  limitation  of  the  radiative  transfer  approach  has  been  identified:  use  of  a 
narrow  bandwidth  approximation.  This  assumption  is  expected  due  to  the  use  of  a 
constant  backscatter  coefficient  with  respect  to  frequency.  However,  using  the  full,  two- 
frequency  mutual  coherence  function,  it  may  be  possible  that  the  impulse  response 
approach  can  be  extended  to  broader  bandwidth  pulses  in  addition  to  broader  beamwidth 
antenna  patterns.  This  premise  is  still  under  investigation. 

As  a  simple  example,  a  simulation  for  a  layered  media  with  a  constant  propagation 
speed  and  constant  extinction  coefficient  in  each  layer  was  performed.  The  results  are 
shown  in  the  following  figures.  The  assumed  parameters  of  the  radar  system  are  as 
follows: 

•  Waveform:  Square  Pulse  with  5  ns  pulse  length 

•  Antenna:  Gaussian  pattern,  1  to  10  degree  beamwidth 

•  100,000m  range  to  the  surface  at  nadir  pointing 
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The  media  is  assumed  to  have  the  following  bulk  propagation  properties  (arbitrary 
estimates  since  no  data  was  available)  and  Gaussian  statistics: 

1 .  The  foliage  layer 

•  20  meter  thickness 

•  effective  extinction  coefficient  as  noted 

•  group  velocity  as  noted 

•  variance  of  the  heights:  0.1  m 

•  backscatter  to  forward  scatter  cross-section,  Cf  /<Jb :  500 

•  absorption  to  total  scattering  cross-section,  aa  /at :  0.9 

2.  The  trunk  layer 

•  5  meter  thickness 

•  effective  extinction  coefficient  as  noted 

•  group  velocity  as  noted 

•  backscatter  to  forward  scatter  cross-section,  <Tf  /g\,  :  500 

•  absorption  to  total  scattering  cross-section,  cra  /cjt :  0.9 

3.  The  ground  layer 

•  Perfect  electric  conductor  (PEC) 

•  variance  of  the  heights:  0. 1  m2 

Figure  12  shows  the  simulated  returned  waveform  for  three  different  antenna 
beamwidths:  1,  5,  and  10  degrees.  As  an  aide  to  the  understanding  of  these  waveforms, 
the  components  of  the  composite  wave  for  the  one-degree  case  are  shown  in  Figure  13; 
here,  the  composite  waveform  is  the  simple  superposition  of  the  volume  and  surface 
scattered  waveforms.  As  you  can  see  from  Figure  13  for  the  one-degree  beamwidth,  the 
surface  return  begins  after  the  volume  return  has  decreased;  hence  in  the  composite 
waveform,  the  surface  and  volume  return  are  still  somewhat  separable.  The  early  return 
of  the  composite  waveform  is  due  to  the  volume  response  and  its  tail  is  attributable  to  the 
surface  return.  However,  from  Figure  12,  for  the  five-degree  beamwidth  case,  this 
distinction  is  less  noticeable,  whereas  for  the  ten-degree  beamwidth  case,  the  surface  and 
volume  returns  become  indistinguishable. 
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Similarly,  in  Figure  14  and  Figure  15,  we  see  waveforms  that  are  derived  using  the 
same  parameters  with  the  exception  of  propagation  speed  in  the  foliage  and  trunk 
regions:  these  have  been  reversed.  In  this  case,  since  the  electromagnetic  depth  of  the 
media,  commonly  referred  to  as  the  optical  depth,  is  much  shorter,  the  surface  and 
volume  returns  are  easily  separable.  Again,  Figure  15  explicitly  shows  the  components  of 
the  composite  waveform  for  a  one-degree  antenna  beamwidth. 


46 


Normalized  Magnitude 


Normalized  Waveform 
(Beamwidth  Variation) 


Beamwidth 

- 1  deg 

- 5  deg 

- 10  deg 


time  (s) 


Figure  12:  Composite  waveform  for  1, 5  and  10  degree  antenna  beamwidths 
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Figure  13:  Components  of  the  1-degree  case:  Volume  +  Surface  Returns 
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Figure  14:  Composite  waveform  for  1, 5  and  10  degree  antenna  beamwidths 
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Figure  15:  Components  of  the  1-degree  case:  Volume  +  Surface  Returns 
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4.2  Modified  Single  Scatter  Theory  for  a  Single  Scatterer  over  a 
Rough  Surface 

4.2.1  Introduction 

As  previously  discussed  the  radiative  transfer  approach  has  simulated  the  foliage- 
surface-foliage  scattering  event  with  only  a  first  order  interaction  in  addition  to  some 
limiting  assumptions.  The  review  of  the  first  order  multiple  scattering  theory  in  Appendix 
A  has  shown  not  only  the  equivalence  of  this  method  to  the  radiative  transfer  method 
when  the  same  assumptions  are  employed,  but  also  a  way  to  generalize  these  results  to 
broader  bandwidth,  broader  beam  systems  and  less  restrictive  scattering  functions. 
Consequently,  in  this  section  we  hope  to  explore  the  extension  of  first  order  multiple 
scattering  to  foliage  over  a  rough  terrain.  This  extension  will  not  only  help  remove  some 
of  the  restrictions,  but  also  provide  a  method  to  examine  the  validity  of  the  first  order 
interaction  with  the  surface.  This  comparison  is  possible  due  to  two  reasons.  First,  using 
the  first  order  multiple  scattering  results,  we  can  produce  results  with  multiple 
interactions  between  the  foliage  and  the  surface.  Second,  for  a  single  scatterer  over  a 
rough  surface,  we  can  directly  compare  these  results  with  exact  numerical  results  as 
presented  in  the  following  section.  The  work  in  this  section  is  still  under  development. 

Appendix  A  has  presented  the  first  order  multiple  scattering  result  for  scattering  from 
a  volume.  This  is  based  on  single  scatter  theory  and  is  equivalent  to  the  radiative  transfer 
approach  previously  outlined.  A  second  important  scattering  event,  in  addition  to  the 
backscatter  due  to  the  volume  alone,  is  the  return  from  the  surface.  There  are  three 
possible  mechanisms  that  may  be  responsible  for  this  process,  see  Figure  16 

1 .  Direct  return  from  the  surface  in  which  no  interaction  with  the  volume  scatterers 
has  occurred. 

2.  A  scattering  event  involving  first  a  member  of  the  volume  of  scatterers  and  then 
interaction  with  the  surface 

3.  A  scattering  event  involving  first  scattering  from  the  surface  and  then  interaction 
with  a  member  of  the  volume  of  scatterers. 
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The  first  case  corresponds  to  free  space  scattering  by  a  rough  surface  that  was  discussed 
in  Appendix  A,  or  the  radiative  transfer  result  discussed  in  Chapter  3.  The  second  and 
third  cases  have  not  been  discussed.  However,  since  the  first  case  has  already  been 
discussed  and  the  second  and  third  cases  are  equivalent  by  reciprocity,  we  need  only 
examine  the  second  case:  the  scatter  from  the  vegetation  to  the  surface  and  back  to  the 
radar.  If  the  scattering  mechanism  for  foliage-to-surface-to-foliage  is  considered 
significant,  this  will  be  a  simple  extension  of  the  results  that  follow. 


Figure  16:  Particle  and  Surface  Scattering  Events 
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4.2.2  Development  of  the  Model 

Given  the  scattering  pattern  of  the  scatterer  and  the  two-frequency  mutual  coherence 
function  of  the  scatterer,  the  power  density  scattered  toward  the  surface  can  be  quantified 
as  follows,  see  Appendix  A  [Ishimaru,  1997]: 


IMsA)  =  111  I  H^’t-fje^ 


volume  of  L“°° 
scatterers 


IT-*-* 


ej(°2td©. 


pdV  (1) 


where  rs  is  the  distance  from  the  antenna  to  the  scatterer  and  the  product  of  the  transfer 
function  and  the  complex  envelope,  Uj(co),  of  the  incident  waveform  is  given  by 


A(t,con)  =  U^coJFCo^k^kOgCo);^)^ - expj-  rsp(tft)dR}  (2) 

rs  ' 

and  the  unit  vectors,  ks  =  (0S, <()s) ,  k;  =  (0;,^),  represent  the  incident  and  scattered 
directions  with  respect  to  the  scatterer,  respectively.  The  function  g(ca,k;)  is  related  to 
the  antenna  pattern;  F(con;ks,ki)is  the  Fourier  Transform  of  the  scattering  function 
f(ks,k;),  and  Uj((Dn)  of  the  complex  envelope  of  the  incident  pulse  at  a  given 
frequency.  See  Appendix  A  for  additional,  more  detailed  definitions. 

The  expression  (1)  describes  the  pulse  shape,  amplitude  and  the  scattering  pattern  of 
the  scattered  power  waveform  due  to  a  volume  of  scatterers,  each  at  a  variable  distance  rs 
from  the  antenna.  The  observation  location  is  on  the  terrain  at  a  distance  rg  from  each 
scatterer.  See  Figure  17.  If  the  narrow  band  approximation  can  be  used,  this  expression 
can  be  simplified  to 


-  ffl  MM.>  Pt[«-7 

volume  of  (471 )  Ts  T  \  C  J 

scatterers  * 
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exp  {-  Psp{crt>dR  }exp {-  p(at)dR  }  pdV 


51 


where abi(ks,kj)  is  the  bistatic  scattering  cross  section  of  the  scatterer  assuming  a 

scattered  waveform  in  the  ks  direction  and  given  an  incident  field  in  the  k;  direction  at  a 
range  rg  with  respect  to  the  scatterer’s  local  coordinates.  Note  that  this  bistatic  cross 
section  can  be  used  due  to  the  narrow  band  assumption  that  results  in  a  particularly 
simple  form  for  the  two-frequency  mutual  coherence  function.  In  general,  this  term  is  the 
integration  of  the  product  of  the  two-frequency  mutual  coherence  function,  the  antenna 
gain  as  a  function  of  frequency  and  the  complex  amplitude  with  its  complex  conjugate 
over  all  frequencies,  see  equation  (2)  and  Appendix  A. 


Figure  17:  Geometry  for  Foliage-Surface-Foliage  Scattering  Event 

Now  we  shall  consider  the  scattered  power  density  from  a  single  scatterer  over  a 
rough  surface.  Using  the  geometry  found  in  Figure  17,  we  can  derive  the  first  order 
multiple  scattering  result  for  this  scatterer.  In  either  case  broadband  or  narrowband,  the 
scattered  power  waveform  created  by  the  scatterer  due  to  its  incident  field  from  the 
antenna  can  be  thought  of  as  a  secondary  source  with  respect  to  the  surface  and  the 
original  source,  the  antenna.  Now  the  surface  has  two  incident  power  waveforms:  the  free 
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space  power  waveform  due  to  the  antenna  and  the  power  waveform  created  by  the 
scatterer.  The  scattered  waveform  from  the  scatterer,  is  described  by  expressions  (1)  and 
(3).  Considering  the  power  waveform  scattered  by  the  scatterer  as  a  second  source, 
l(t;0s,<j>s)  represents  the  product  of  the  transmitted  power  waveform  due  to  the  scatterer, 
PT,s(t),  as  modified  by  an  assumed  “antenna  pattern”  for  the  scatterer,  Gs(0s,(j)s)  .  First 
we  will  re-express  equations  (1)  and  (3)  for  the  specific  case  of  a  single  scatterer. 


Figure  18:  the  geometry  of  the  scatterer 

If  we  consider  a  single  scatterer  under  the  narrow  band  approximation,  the 
expression  for  the  power  density  incident  on  the  scatterer  is  found  under  the  first  order 
multiple  scattering  theory  to  be 

■  A )  exp{_  jfe  }  pit  -if 

47rrs  V  c 

here,  PT(t)  is  the  radar’s  transmitted  waveform  that  has  been  attenuated  by  the 
exponential  factor  and  weighted  by  the  radar’s  antenna  gain  Ga(0a,<j>a)  in  the  direction 

of  the  scatterer,  (0a ,  (j)a ) .  These  angles  describe  the  angle  of  the  scatterer  with  respect  to 
the  antenna  in  the  antenna’s  reference  frame  and  rs  is  the  distance  from  the  antenna  to  the 
scatterer.  Employing  the  bistatic  radar  cross-section,  crbj(ks,kj),  the  waveform  scattered 
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A  A 

in  the  ks  or  (0s,<j>s)  direction  given  an  incident  field  in  the  kj  direction  with  respect  to 
the  scatterer’s  local  coordinates  can  be  written  as 


UWM  =  {G,f‘T-;)  «p{-  J'‘k.(n)dM }  PTft  - 5.1}  CTbi(k„k,) 


The  “scattering  pattern”  of  the  scatterer  in  the  direction  from  the  scatterer  to  the  terrain, 
(0g,  (|)g),  is  the  analogous  normalized  antenna  pattern  of  the  new  source  with  a  normalized 
antenna  pattern  given  by 


Gs(Mg)  = 


MMi) 

_max 

abi 


(5) 


and  the  analogous  transmitted  power  density  from  the  scatterer  is  given  by 

PT.s(t)=  e*p{-  I'K  (n)dM  }  ar  P,  f  t  -  -1  (6) 

47trs  V  c  J 

Note  that  the  product  of  (5)  and  (6)  is  equal  to  (4).  Recall  that  the  angles  using  subscript 
“a”  refer  to  the  angle  between  the  antenna  and  the  scatterer  in  the  antenna’s  local 
coordinate  system  and  rs,  refers  to  the  slant  range  from  the  antenna  to  the  scatterer;  see 
Figure  17.  In  the  general  broadband  case,  these  quantities  are  not  so  easily  separated 
since  they  remain  buried  in  the  inverse  transforms.  This  matter  requires  further 
investigation. 

Next,  we  consider  the  foliage-surface  scattered  power  collected  by  the  radar.  From 
the  discussions  involving  the  flat  surface  impulse  response,  the  solution  for  the  average 
incoherent  response  of  a  scatterer  over  a  rough  surface  can  be  developed.  We  begin  with 
an  expression  for  the  returned  power  waveform  which  is  similar  to  that  given  in  Chapter 
2  describing  the  impulse  response  method  for  a  rough  surface  in  free  space, 


PR(t)  = 


f  2(ro-^(x,y)sec0)" 

12  «2xrT  1 

±  f  f_V _ _ )_ 
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In  this  case,  the  transmitter  and  receiver  are  at  the  same  location  and  have  the  same 
antenna  pattern;  the  slant  range  from  antenna  to  surface  is  simply  traversed  twice.  When 
the  scatterer  acts  as  a  source,  however,  and  the  radar  antenna  collects  energy,  the 
expression  must  change  to  reflect 


1.  the  different  distances:  from  the  scatterer  to  the  surface  element,  rg  -  ^(x,  y)  sec  0S 
and  from  the  radar  antenna  to  the  surface  element,  r0  -£(x,y)sec0  ;  these  range 
variations  affect  the  delay  time  as  well. 

2.  the  different  “antenna  patterns”:  from  the  scatterer,  Gs(0s,<j>s)  with  respect  to  the 
surface  and  from  the  radar  antenna,  Ga  (0,  (j>) ,  with  respect  to  the  surface,  each  in 
its  own  coordinate  system. 

/V 

3.  The  bistatic  radar  cross  section  per  unit  area  with  an  incident  angle,  k;  or 

A 

(0i , (j); )  and  scattered  angle  ks  or  (0s,<j)s)  with  respect  to  the  antenna’s  coordinate 
system,  see  Figure  18. 

Hence,  the  expression  for  the  incoherent  power  returned  from  the  surface  as  a  function  of 
time  due  to  the  scattered  energy  in  the  foliage-surface  interaction  can  be  written  as 


1  ^  co  2  7C 
•  *  re 


^T,S  t 


(rg-5(x,y)sec0s)  (r0  -  foe,  y)  sec  0) 
co  co 


P"(tl  (4*)1  !  I  (r  -  tfx, y) sec  9S  f  (r0  - 5(x,  y)  sec  e)! 


•Gs(e„+,)G.(e,*)  Pd*dP 


note  that  the  incident  power  waveform  is  that  from  the  scatterer,  Pt,s  defined  by  equation 
(6).  Simplifying  this  result  under  the  assumption  that  the  surface  height  is  negligible  in 
the  amplitude  terms,  and  recognizing  the  equivalence  of  the  radar  coordinates  and  the 
polar  coordinates  defined  on  the  surface,  we  find 
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The  integration  over  the  surface  is  performed  in  the  antenna  coordinates.  Assuming 
the  angles  above  can  be  measured  from  the  mean  plane,  a  relationship  between  the 
“scatterer-to-ground”  angles,  subscripted  “s”,  and  the  “antenna-to-ground”  angles,  not 
subscripted,  can  be  found  with  respect  to  the  antenna  frame  of  reference  in  the  following 
symbolic  manner.  Let 


0S  =  fs(0)  =  tan’ 


<f>s  =  gs(<t>)  =  ±tan 


-if  (x-x0)2+(y-y0)2 


x~xo 
v  y-y0 


In  addition  the  slant  range  from  the  scatterer  to  the  terrain  can  be  found  in  terms  of  the 
antenna  coordinates.  Given  that  the  scatterer  is  located  at  the  coordinates  (x0,  yo,  hs),  the 
distance  from  the  scatterer  to  the  surface  is  expressed  as  follows,  see  Figure  1 7 

rg  =  (x-x0)2+(y-y0)2+  hs2 

with  (x,  y)  =  ( r0  sin  0  cos  <f> ,  r0  sin  0  sin  <|>) 


Hence,  establishing  a  reference  frame  on  the  antenna,  we  express  all  angles  relative  to  the 
spherical  coordinates  centered  at  the  antenna.  Consequently,  expressed  entirely  in  the 
antenna’s  frame  of  reference,  the  power  return  waveform  becomes 
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where  each  angle  and  distance  is  now  measured  in  terms  of  the  antenna/radar  coordinates. 
Averaging  these  results  by  the  surface  heights,  the  average  incoherent  power  returned 
becomes 
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where  the  probability  density  function  (pdf)  for  the  surface  heights  has  been  given  by 
p.(^) .  With  a  change  of  variables  for  the  heights  given  by 


|(x  y)  =  £(x,y)(  sec(f5(e))+  sec  6) 
cn 


and  the  corresponding  transformation  of  the  pdf  given  by 
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the  average  returned  power  can  be  re-written  as 
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the  returned  average  incoherent  power  can  be  expressed  in  a  convolutional  form,  with  the 
symbol  <£>  denoting  convolution. 
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Taking  advantage  of  the  properties  of  the  delta  function, 
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Next,  we  invoke  the  narrow  beam  approximation  for  the  dependence  of  the  argument  for 
the  pdf  and  we  approximate  0  with  the  boresight  angle,  0q. 


/,  \ _  c0  Cq  t 

P|  sec(fs  (0O ))  +  sec  0O  P \  sec(fs  (0O ))  +  sec  0O 


Finally,  the  expression  for  the  average  power  returned  from  the  scattering  event:  antenna- 
to-scatterer-to-ground-to-antenna,  can  be  represented  by  the  convolutional  product 

(Pr («)>  3PT.s(l)®P5(t)®PFs.s(t)  (12) 


where 


Pfs.s  (1) 


s  ft,  (fr(r0)  +  r0)' 

y2  OO  2lt  i  p 

-  (4 nf  II  (sec(fs(e))+sec6) 


•  Gs(f,  t0).  g,  W)G,  (0,  <t»)cr0  (f,  (0),  gs(4>>6.4>)  r„  d+dr. 


Note  that  the  height  of  the  single  scatterer,  hs,  and  its  orientation,  Q,  can  be  made  random 
variables  as  well;  this  formulation  has  yet  to  be  investigated.  The  random  orientation  of 


the  scatterer  will  require  that  the  boresight  angle  of  Gs,  the  scatterer’s  “gain  pattern,”  be  a 
random  variable.  A  scatterer  with  a  random  height  will  result  in  the  distance  rg  and  the 
associated  angles  to  become  random  variables.  The  final  average  return  power  from  the 
surface  due  to  a  single  scatterer  will  then  require  that  we  average  only  the  modified  flat 
surface  impulse  response,  Pfs,s(0-  Hence  the  averaged  incoherent  power  returned  will  be 
of  the  form 
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where  hs  =  a  random  height  for  scatterer  N 

Q  =  a  random  orientation  for  scatterer  N 


and(«)a  is  interpreted  as  ensemble  averaging  with  respect  to  the  random  variable  a 

Now  the  return  from  the  surface  and  the  scatterer  combination  is  a  superposition  of 
the  return  from  the  volume  and  that  from  the  surface.  Consequently,  the  incoherent 
response  from  N,  randomly  oriented,  scatterers  at  random  heights  can  be  expressed  as 

(P(t))=  P„(t)®P5(t)®PK(t)  +£PT.s(t)®P;(t)®(PFS.s(t))  (14) 


The  first  term  is  the  incoherent  power  received  directly  from  the  scatterers  as  described  in 
the  radiative  transfer  section  or  the  first  order  multiple  scattering  description.  The  second 
term  is  the  response  for  N,  randomly  oriented  scatterers  at  random  heights  above  a 
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random  surface.  A  third  term  may  be  added,  if  it  is  expected  that  there  will  be  significant 
return  directly  from  the  surface  without  interaction  from  the  scatterers. 

4.2.3  Conclusions  and  Future  Efforts 

In  this  section,  an  outline  for  predicting  the  average  incoherent  power  returned  from  a 
single  scatterer  and  N  randomly  oriented  scatterers  at  a  random  height  above  a  rough 
surface  has  been  presented.  Implemented,  this  model  will  accomplish  the  following  with 
respect  to  the  radiative  transfer  model. 

1.  When  N  scatterers  are  used,  the  radiative  transfer  results  should  be  duplicated  if  a 
single  interaction  is  considered  along  with  the  narrow  beam,  narrow  bandwidth  and 
restricted  scattering  pattern  assumptions  as  detailed  in  Chapter  3,  the  radiative 
transfer  approach. 

2.  The  above  restrictions  may  be  relaxed  in  the  first  order  multiple  scattering  approach 
in  order  to  judge  their  effect;  thus  supporting  or  further  generalizing  the  radiative 
transfer  result  of  Chapter  3. 

•  The  scattering  pattern  may  be  chosen  more  realistically 

•  The  beamwidth  of  the  illuminating  antenna  and  the  bandwidth  of  the  signal  may  be 
increased  (with  the  possible  sacrifice  of  the  convolutional  form) 

•  Including  the  higher  order  interactions  between  the  foliage  and  the  surface  may  check 
the  single  passage  assumptions. 

•  If  the  radiative  transfer  results  prove  to  be  too  restrictive  through  the  implementation 
of  these  generalizations,  this  model  will  serve  as  the  next  possible  approach. 

3.  Simulation  of  a  single  scatterer  above  a  rough  surface  may  be  checked  via  exact 
numerical  calculations  as  outlined  in  the  next  section  of  this  chapter,  thus  checking 
the  first  order  multiple  scattering  approach  (with  higher  order  interactions  between 
the  surface  and  the  scatterer)  and  ultimately  the  radiative  transfer  model. 
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The  implementation  of  a  single  scatterer  placed  deterministically  above  a  rough 
surface  must  be  accomplished  first.  These  results  will  be  compared  with  the  exact  results 
as  discussed  in  the  next  section.  The  formal  implementation  of  a  randomly  oriented  and 
randomly  positioned  particle  must  then  be  formulated.  Averaging  the  new  Flat  Surface 
Impulse  Response  will  be  attempted.  Once  this  accomplished,  we  extend  this  technique 
to  N  scatterers  as  discusses  in  the  previous  section.  Using  the  Foldy-Lax-Twersky 
integral  equation,  we  can  establish  the  effective  media  for  the  first  order  multiple 
scattering  results  and  obtain  a  result  that  is  similar  to  the  Distorted  Wave  Bom 
Approximation  (DWBA);  consequently,  these  results  may  be  compared  with  those 
obtained  by  Lang  [1981].  Finally,  the  broadband  approach  must  be  fully  explored.  The 
challenges  in  this  extension  include  keeping  our  results  in  the  fully  convolutional  form 
obtained  for  the  narrowband  case.  Most  likely,  we  will  choose  a  scattering  function  that  is 
simple  and  practical  such  as  that  presented  by  Schwering  [1985] 
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4.3  Exact  Solution  to  the  Single  Scatterer  above  a  Rough  Surface 

The  impulse  response  model,  like  most  radiative  transfer  models,  does  not  account  for 
any  interaction  between  the  scatterers  (foliage)  and  the  boundary  (surface).  In 
establishing  a  range  of  validity  for  this  assumption,  a  measure  and  threshold  of  “no 
interaction”  must  be  established.  Once  this  measure  is  established,  numerous  simulations 
of  the  exact  scattered  field  must  be  examined  in  order  to  verify  this  assumption  over  a 
large  parameter  space  including 

•  the  scatterer’ s  size  normalized  to  wavelength 

•  the  scatterer’s  separation  from  the  rough  surface  normalized  to  wavelength 

•  the  scatterer’s  orientation  (if  it  is  not  circular) 

Consequently,  we  must  assess  the  magnitude  of  the  contribution  of  multiple  scattering 
interactions  between  the  scatterer  and  the  rough  surface.  One  approach  to  establishing  the 
measure  of  significant  interaction  would  be  to  include  each  level  of  multiple  scattering 
and  measure  its  contribution  to  the  exact  solution.  Hence,  we  begin  with  the  assumption 
that  the  surface  and  the  scatterer  do  not  interact.  Next,  we  assess  the  correction  for  a 
single  scatter  interaction. 

In  order  to  verify  this  assumption  of  independent  scattering,  an  “exact”  numerical 
model  has  been  created  using  the  method  of  moments  (MOM);  the  numerical  solution  for 
the  currents  on  the  scatterer  and  the  rough  surface  accounts  for  all  orders  of  interactions. 
After  the  problem  is  cast  into  the  proper  integral  equation,  the  geometry  is  discretized  in 
preparation  for  a  solution  via  the  familiar  Method  of  Moments  (MOM).  Specifically,  the 
Method  of  Ordered  Multiple  Interactions  (MOMI)  has  been  modified  and  is  implemented 
as  a  solution  method.  This  will  be  described  in  Section  14.34.3.1.  Once  the  currents  are 
found  from  the  integral  equation,  the  scattered  fields  can  be  simply  found  using  the 
proper  radiation  integral;  the  far-field  formulation  has  been  used.  A  brief  description  of 
the  MOMI  as  originally  applied  to  rough  surfaces  can  be  found  in  Chapter  2. 
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4.3.1  Extension  of  the  MOMI  to  Closed  Bodies 


The  MOMI  method,  as  discussed  in  Chapter  2,  is  a  solution  method  for  the  MOM 
derived  matrix  equation  of  the  following  form: 

\j/  =  \|/“+P\|/  (1) 

where  P  is  a  propagator  matrix,  vy  is  an  unknown  scalar  field  and  vj/inc  is  the  known 
incident  field.  In  developing  MOMI  to  analyze  scattering  from  extended  rough  surfaces, 
the  self-interaction  terms  Pjj  were  neglected  [Kapp  and  Brown,  1996].  The  propagator 
matrix  (P)  was  thus  decomposed  into  lower  triangular  (L)  and  upper  triangular  (U) 
matrices,  each  having  zero  entries  along  the  diagonal, 

P  L  +  U  (2) 

After  a  few  simple  manipulations,  this  decomposition  led  to  the  MOMI  matrix  equation 
given  in  Chapter  2.  However,  consistent  discretization  of  (1)  requires  that  the  diagonal 
elements  Pjj  be  retained  [Toporkov,  1998].  This  modification  can  be  incorporated  in 
several  ways.  In  applying  MOMI  to  integral  equations  having  singular  kernels,  it  has 
been  found  that  optimal  convergence  properties  are  obtained  by  decomposing  the 
propagator  matrix  as 

P  — »  L  +  D  +  U  (3) 

where  D  is  a  diagonal  matrix  with  D  =  Pjj.  Physically,  maintaining  the  self  interaction 
terms  in  D  separate  from  (L)  and  (U)  provides  better  convergence  properties  when 
applying  the  method  to  integral  equations  having  singular  kernels  because  these 
equations  exhibit  strong  coupling  between  oppositely  directed  fields  on  the  surface  of  a 
scatterer  [Adams  and  Brown,  1997], 

The  decomposition  (3)  leads  to  the  matrix  equation 

V=(D-U)-,D(D-L)-‘lf“  +  P„V  (4) 

where  D  =  I  -  D  and  the  MOMI  propagator,  Pm,  is  defined  as 

PM  =  (D  -  U)"1  D  (D  -  L)-1  LD"1  U 
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(5) 


Neumann  iteration  of  (4)  yields  the  candidate  solution 


>|/  =  XV(D-Ur'D(D-L)-V“  (6) 

11=0 

which  is  the  same  as  equation  (28)  of  [Kapp  and  Brown,  1996]  under  the  substitution 
D->  I .  Since  D  =  I  +  0(Ax),  the  convergence  properties  of  (6)  are  essentially  unchanged 
from  those  of  equation  (28)  of  [Kapp  and  Brown,  1996].  When  convergence  occurs,  the 
candidate  solution  (6)  converges  to  the  exact  solution  of  (1).  As  discussed  below,  even  in 
cases  for  which  the  infinite  series  does  not  converge,  the  first  few  terms  of  (6)  can  still 
provide  a  good  approximation  to  the  actual  solution. 

The  MOMI  series  (6)  provides  a  very  robust  and  rapidly  convergent  solution  to  the 
MFIE  for  scattering  from  extended  rough  surfaces  in  two  dimensions.  The  series  has 
never  been  observed  to  diverge.  These  desirable  properties  have  been  attributed  to  the 
manner  in  which  the  MOMI  series  re-sums  the  multiple  scattering  terms  present  in  the 
Neumann  series  for  the  original  integral  equation  (1).  The  Bom  term  in  the  MOMI  series, 
(D-U)"ID(D-L)'1v|/m%  includes  the  contributions  to  the  current  due  to  all  orders  of 
continuous  forward  scattering  (D  -  L)'1 ,  all  orders  of  backscattering,  (D  -  U)"1 ,  and  one 
order  of  interaction  between  the  backward  and  forward  traveling  waves  on  the  surface 
(resulting  from  the  multiplication  of  these  operators).  Thus,  the  Bom  term  includes 
interactions  of  up  to  order  N.  The  largest  effect  neglected  by  the  order  zero  iterate  of  the 
MOMI  series  is  that  of  a  wave  which  twice  changes  directions  on  the  rough  surface 
before  again  interacting  with  the  currents  on  the  surface  -  a  triple  scattering  event. 

For  this  reason,  the  ordering  of  the  unknowns  in  the  original  matrix  equation  (1) 
can  have  a  drastic  effect  on  the  convergence  of  the  MOMI  series.  This  is  in  contrast  to 
the  Neumann  series  for  (1)  whose  convergence  properties  are  independent  of  the  manner 
in  which  the  unknowns  are  ordered  in  the  matrix  equation.  A  different  ordering  of 
unknowns  in  the  MOMI  series  will  result  in  the  summation  of  different  multiple 
scattering  terms.  In  the  case  of  a  random  ordering  of  the  unknowns  in  the  original  matrix 
equation  (1)  for  the  rough  surface  scattering  problem,  the  number  of  MOMI  iterations 
required  to  converge  to  a  given  error  tolerance  can  be  orders  of  magnitude  larger  than  in 

64 


the  case  of  the  physically  based  forward-backward  ordering.  It  is  not  immediately  clear 
how  the  unknowns  in  (1)  should  be  ordered  for  the  application  of  MOMI  to  closed  body 
scattering  problems.  In  the  case  of  elliptical  cylinders,  at  least  two  ordering  schemes 
incorporate  important  physical  aspects  of  the  scattering  problem.  These  methods  of 
ordering  the  unknowns  in  the  matrix  equation  are  illustrated  in  Figure  19.  An  ordering 
which  is  sequential-in-<()  (SIP)  produces  an  iterative  series  that  mimics  the  progression  of 
creeping  waves  around  the  surface  of  the  cylinder.  An  alternative  approach  is  one  that  is 
sequential-in-x  (SIX).  This  ordering  results  in  a  MOMI  series  for  the  closed  body 
problem  which  is  somewhat  analogous  to  the  forward-backward  approach  used  in  [Kapp 
and  Brown,  1996]. 


SIP  Ordering  SIX  Ordering 

Figure  19  :  Ordering  of  the  Unknowns 
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4.3.2  A  Combined  Field  Formulation 


Given  this  understanding,  we  next  consider  a  formulation  of  the  scattering  problem 
that  does  not  give  rise  to  a  singular  or  nearly  singular  integral  equation.  In  the  following 
we  consider  a  combined  field  integral  equation  (CFIE)  representation  [Mautz,  1978].  The 
CFIE  is  a  linear  combination  of  the  MFIE  and  the  EFIE  as  indicated  below. 

aEFIE  +  MFIE  =CFEEa 

While  the  CFIE  can  be  used  to  provide  a  unique  solution  to  the  scattering  problem, 
the  use  of  MOMI  as  formulated  with  a  combined  field  description  of  the  scattering 
problem  introduces  additional  difficulties  associated  with  the  kernels  of  the  EFIEs.  The 
EFIE  kernel  for  the  TE  problem  is  simply  the  Green's  function  of  the  Helmholtz  equation. 
For  TM  scattering  the  kernel  function  is  the  second  normal  derivative  of  the  Green's 
function.  In  applying  the  MOMI  series  to  scattering  from  dielectric  surfaces  it  has  been 
found  that  the  singularities  of  these  kernel  functions  produce  strong  coupling  between 
oppositely  directed  fields  [Adams  and  Brown,  1997].  The  singularity  present  in  the  EFIE 
for  TM  scattering  is  particularly  strong  and  requires  that  a  modified  form  of  (4.3. 1-6)  be 
used.  The  singularity  of  the  EFIE  kernel  for  TE  scattering  is  much  weaker  and  is 
therefore  more  amenable  to  the  MOMI  series  solution  technique.  For  this  reason,  in  the 
following  we  investigate  the  application  of  MOMI  to  the  CFIE  for  TE  scattering  only. 


4.3.3  Selecting  an  optimal  CFIE:  A  multiple  scattering  approach 

The  electric  field  integral  equation  (EFIE)  for  TE  scattering  from  a  PEC  object  is 

0  =  £inc-  ff—  GdS0  (1) 

sJano 


The  CFIE  for  the  TE  case  is  obtained  by  adding  this  to  the  MFIE  for  this  problem  using 
the  complex  constant  a.  This  leads  to 


dE 

dn 


-2aEmc  +  2 


aEmc 

an 


(2) 
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where 


Ka- 


f  r  5Gn 
aG  +  — 
dn  ) 


(3) 


Discretization  of  this  equation  is  discussed  in  [Adams  and  Brown,  1998b].  The  resulting 
matrix  equation  can  be  put  in  the  form  of  (4.3. 1-1).  The  corresponding  MOMI  series  is  of 
the  same  form  as  (4.3. 1-6). 

The  CFIE  (2)  is  guaranteed  to  have  a  unique  solution  whenever  a  is  complex.  This 
requirement  provides  significant  freedom  in  the  choice  of  a.  We  further  constrain  a  by 
requiring  that  it  provide  optimal  convergence  properties  for  an  arbitrary  incident  field. 
This  corresponds  to  the  value  of  a  which  minimizes  the  maximum  modulus  eigenvalue  of 
Pm- 

A  physically  intuitive  way  of  determining  this  optimal  choice  is  to  minimize  the 
contributions  to  the  total  field  in  (2)  which  are  due  to  the  integral  term.  Observe  that  if  a 
is  chosen  such  that  Ka  =  0  then  no  iteration  is  required  to  obtain  the  exact  solution  to  the 
scattering  problem.  This  is  physically  interpreted  as  the  case  of  zero  multiple  scattering, 
an  example  of  which  occurs  in  the  MFIE  (a  =  0)  formulation  of  scattering  from  an 
infinite  PEC  planar  surface.  In  this  case 


dn 


and  the  integral  term  in  (2)  provides  no  contribution  to  the  total  surface  current.  This 
results  for  the  flat  surface  scattering  problem  because  the  magnetic  field  radiated  to  an 
observation  point  on  the  surface  by  sources  which  are  also  on  the  surface  has  no 
tangential  component.  Thus,  the  MFIE  for  which  MOMI  was  originally  developed  can  be 
seen  as  the  specialization  of  (2)  to  the  rough  surface  case  where  a  is  selected  such  that  K« 
=  0  in  the  unperturbed  geometry  of  the  rough  surface  scattering  problem. 

If  possible,  the  optimal  choice  for  a  would  in  general  be 


J_5G 
G  dn 


(4) 
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as  this  would  give  K«  =  0  in  all  cases.  Because  this  choice  is  not  generally  possible,  we 
instead  consider  asymptotically  determined  estimates  of  a.  In  addition  to  providing 
optimal  integral  formulations,  these  asymptotic  estimates  provide  information  on  how  the 
optimal  value  of  a  depends  on  the  size  and  shape  of  the  scatterer.  To  simplify  the 
following  analysis,  we  consider  only  the  case  of  circular  cylinders. 

For  small  cylinder  radii  (ka  «  1)  the  normal  derivative  of  the  electric  field  on  the 
surface  is  approximately  constant.  Thus 

rrd 8E 


and  the  contribution  of  the  integral  term  to  the  total  field  in  (2)  is  minimized  by  choosing 


a  =  - 


(dG/dn) 


where  (•)  denotes  an  average  of  the  source  point  over  the  surface  of  the  cylinder.  Using 
the  small  argument  forms  of  the  zero  and  first  order  Hankel  functions 


H"»(z)~1-£  lnffl+y 
n  12 


where  y  =  0.5772156649 ,  the  Euler-Masheroni  Constant 


substituting  above  gives 


H[2)  (z)  ~ 


j2  /  f  j2y>)  j2  [  .  f  <f>-.<|>0 

a  In  kasm  - — — 

a\l  n  J  n  [  l 2 


■  —  [j7r  +  21n(ka)-  ln2  +  y] 


which  is  seen  to  be  inversely  related  to  the  size  of  the  cylinder. 
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For  large  cylinder  radii  (ka  »1)  the  approximation  (5)  is  not  valid.  Instead  we 
attempt  to  minimize  the  contribution  of  the  integral  term  appearing  in  (2)  by  choosing  for 
a  the  average  of  (6),  i.e., 

/dG/dn\ 

\  G  / 


Substituting  the  large  argument  forms  of  the  zero  and  first  order  Hankel  functions 

Hq2)  (z)~  J^e-j(z-n/4) 

V  TCZ 

H{2)  (z)~1/^e-j(z-35t/4) 

V  7TZ 


the  simplified  a  becomes1 

j2k  =  j4 
n  X 


a  ~  jk 


sin 


<1>— 4>0 


(8) 


Unlike  the  small  ka  limit  for  which  the  optimal  choice  of  a  was  found  to  vary  with  a,  the 
value  of  a  suggested  by  (8)  is  independent  of  the  cylinder's  size. 

The  physical  significance  of  the  value  of  a  given  by  (8)  is  understood  by  observing 
that,  on  the  surface  of  the  cylinder  in  the  TE  problem 

3E 

—  y  =  jcopn  x  H  =  jcopHBny 
on 

Multiplying  this  equation  through  by  a'1  and  inserting  the  asymptotic  value  of  a  given  by 
(8)  prior  to  averaging,  we  have 


Although  the  averaging  performed  in  (7)  has  been  performed  over  all  <j>0,  the  contributions  to  the  average 
from  points  near  the  observation  point  (for  which  the  large  argument  approximations  of  the  Hankel 
functions  are  not  valid)  are  of  order  ka'1  when  ka  is  large  and  do  not  significantly  affect  the  result  for  a. 
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where 


1  dE  =  jcop 

a  dn  jk| sin([<j>- 4>0]/2) | 


sin([<t>-<(>0]/2)| 


and  t|  is  the  free  space  impedance.  Recalling  that  a  plane  wave  propagating  in  the  k 
direction  satisfies  the  relation 


r|kx  H  =  -E 


we  see  that  the  choice  of  a  specified  by  (8)  results  from  the  fact  that  the  field  at  an 
observation  point  excited  by  a  distant  source  point  is  locally  planar.  |  sin(  [  cj)  —  <j>0  ]/ 2  )  |  in 
the  denominator  of  r\  (which  is  not  present  above)  arises  because  the  CFIE  imposes  a 
boundary  condition  on  the  tangential  components  of  the  total  field.  For  TE  incidence,  the 
electric  field  is  always  tangential  to  the  cylinder  while  the  component  of  the  magnetic 
field  tangent  to  the  cylinder  varies  with  the  location  of  the  source  point. 

The  optimality  of  the  asymptotic  estimates  for  a  provided  by  (7)  and  (8)  can  be 
evaluated  by  calculating  the  eigenvalues  of  the  resulting  propagator  matrices  Pm-  Figure 
20  and  Figure  21  show  the  magnitude  of  the  largest  eigenvalue  of  the  MOMI  propagator 
Pm  (in  dB)  for  the  TE  CFIE  as  a  function  of  the  complex  constant  a  for  various  cylinder 
radii. 
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Figure  20:  Contours  of  the  largest  eigenvalue  of  Pm  (in  dB)  as  a  function  of  the  complex 

constant  a  for  circular  cylinders 


As  anticipated  by  the  above  asymptotic  analysis,  for  small  radii  the  choice  of  a  which 
yields  the  minimum  maximum  eigenvalue  of  Pm  varies  significantly  with  the  radius  of 
the  cylinder.  The  estimate  of  a  provided  by  (7)  is  good  through  a  =  0.051.  For  the  larger 
radii  cases  of  a  =  11  and  a  =  51  the  value  of  a  which  yields  the  minimum  maximum 
eigenvalue  remains  fairly  constant  at  a  *  j4/l,  which  is  in  good  agreement  with  the 
asymptotic  value  provided  by  (8). 
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Figure  21 :  Contours  of  largest  eigenvalue  of  Pm  (in  dB)  as  a  function  of  the  complex 

constant  a  for  circular  cylinders 

From  these  figures  we  also  notice  that  the  minimum  value  of  the  maximum 
eigenvalue  of  Pm  as  a  function  of  a  increases  as  the  cylinder  radius  increases  from  0.01  A. 
to  5X.  This  occurs  because  for  large  cylinder  radii  we  are  able  to  minimize  the  integral 
term's  contribution  to  (2)  only  by  minimizing  in  an  average  sense.  In  the  small  radii 
limit  the  approximation  in  (5)  allowed  us  to  determine  the  value  of  a  by  minimizing  the 
contribution  from  the  integral  term  itself. 

Figure  22  and  Figure  23  show  the  magnitude  of  the  largest  eigenvalue  of  Pm  for  TE 
scattering  from  elliptical  cylinders  having  axial  ratios  of  a/b  =  8.  For  small  radii  we  see 
that  the  optimal  choice  of  a  changes  significantly  with  the  size  of  the  ellipse.  The 

smallest  maximum  eigenvalue  of  Pm  is  larger  for  small  elliptical  cylinders  than  for  small 
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circular  cylinders  due  to  the  breakdown  of  approximation  (5)  as  the  axial  ratio,  a/b, 
increases. 


Figure  22:  Contours  of  largest  eigenvalues  of  Pm  (in  dB)  as  a  function  of  the  complex 
constant  a  for  elliptical  cylinders  having  a/b  =  8 


Also  note  that  the  optimal  choice  of  a  appears  to  bifurcate  into  two  distinct  regions  in 
the  complex  a-plane  for  small  value  of  a.  This  is  due  to  the  loss  of  rotational  symmetry 
when  (a/b  ^  1)  and  suggests  that  it  may  be  more  appropriate  to  choose  a  as  a  function  of 
position  in  this  case,  i.e.,  a  =  a(p).  For  larger  values  of  a,  the  estimate  in  (8).  derived  for 
circular  cylinders  is  seen  to  provide  an  excellent  choice  in  the  case  of  a/b  =  8. 
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Figure  23:  Contours  of  largest  eigenvalues  of  Pm  (in  dB)  as  a  function  of  the  complex 
constant  a  for  elliptical  cylinders  having  a/b  =  8 


4.3.4  Example  Results  for  TE  Polarization 

Simulations  that  follow  demonstrate  the  investigations  into  an  elliptical  cylinder 
above  a  rough  surface.  A  MOMI  code  was  produced  which  includes  unknowns  on  both 
an  elliptical  cylinder  and  a  rough  surface.  The  coupling  parameter,  a,  on  each  surface 
was  chosen  to  be  the  optimum  for  the  observation  point  on  the  surface.  Consequently,  a 
is  the  asymptotic  value  described  above  on  the  cylinder  and  a  =  0  on  the  surface.  The 
cylinder  was  chosen  to  be  elliptical  in  order  to  resemble  the  foliage  problem.  All  of  the 
simulations  which  follow  use  TE  Polarization  and  A./10  sampling.  Other  parameters 
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depend  upon  the  size  of  the  scatterer  and  its  separation  from  the  surface;  these  include  the 
following 

•  incident  spot  size  =  20a  to  301 

•  surface  length  =  1 001  to  2001 

•  ellipse:  major  axis  =  61  to  101  ,  minor  axis  =  21  to  2.51 

•  height  of  the  ellipse  =  201  to  601  above  the  mean  surface 

These  variables  are  identified  in  Figure  24. 


Figure  24:  Scatterer  over  a  Randomly  Rough  Surface 


The  total  cross  section  of  the  elliptical  cylinder  and  rough  surface  combination  is  plotted 
in  the  figures  that  follow.  The  far  field  form  of  the  Hankel  function  normalizes  this  value: 

&'1; 
v  7tkp 

This  solution  which  includes  the  full  interaction  will  be  referred  to  as  the  “exact 
solution”. 
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In  addition  to  the  total  cross  section  of  the  full  interaction  problem,  the  total  cross 
section  of  various  stages  of  interaction  are  also  included.  First  of  all,  the  simple 
incoherent  addition  of  scatter  power  is  included;  this  case  will  be  referred  to  as 
“incoherent  addition”  (IA)  in  the  examples  that  follow.  This  curve  will  allow  the 
assumptions  of  the  impulse  response  method  (no  interaction)  to  be  compared  with  a  full 
interaction  solution.  Hence,  it  is  equivalent  to  calculating  the  total  cross  of  the  ellipse  and 
the  surface  in  isolation  and  simply  adding  the  resulting  power,  see  Figure  25(a)  and 
Figure  25(b),  respectively.  The  source  in  this  case  is  will  be  referred  to  as  the  “free-space 
incident  field”  and  the  resulting  induced  currents  as  the  “incident  currents.” 


Figure  25:  The  Single  Scatter  Approximation:  (a)  Currents  induced  on  the 
cylinder  in  isolation  (b)  Currents  induced  on  the  surface  in  isolation 

The  next  step  is  the  addition  of  the  first  order  correction  to  both  the  current  on  the 
surface  and  the  ellipse.  Thus,  after  the  current  on  the  ellipse  is  calculated  in  isolation,  its 
radiated  field  is  added  to  the  free-space  field  incident  on  the  surface,  see  Figure  26(a). 
This  results  in  the  simple  correction  of  the  incident  currents  and  the  single  scatter  currents 
due  to  the  ellipse.  In  turn,  this  corrected  current  on  the  surface  is  permitted  to  radiate  and 
induce  a  correction  to  the  current  on  the  ellipse.  This  results  in  a  double  scatter  event,  yet 
is  still  a  first  order  correction  to  the  incident  current  on  the  ellipse,  see  Figure  26(b). 
Finally,  the  composite  system  with  the  first  order  corrected  currents  is  allowed  to  radiate. 
This  result  is  referred  to  as  the  double  scatter  result  in  the  following  example  results. 
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Figure  26:  The  Single  Scatter  Approximation:  (a)  Corrections  to  currents  induced  on  the 
surface  (b)  Correction  to  Currents  induced  on  the  cylinder 


Figure  27:  6  wavelength  ellipse,  60  wavelengths  above  a  Gaussian  rough  surface 
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In  Figure  27  through  Figure  29,  we  see  the  effect  of  separation  on  the  total  cross 
section  for  the  elliptical  scatter/rough  surface  system  by  varying  the  separation  between  a 
6  wavelength  elliptical  cylinder  and  the  surface.  In  these  figures,  the  exact  total  cross 
section  can  be  compared  with  the  incoherent  addition  of  the  surface  and  the  scatterer 
(radiative  transfer  assumption)  and  the  first  order  interaction  between  these  parts.  It  is 
obvious  from  these  figures  that  the  double  scatter  approximation  provides  a  better 
estimate  of  the  total  cross  section  of  the  composite  system  than  the  simple  incoherent 
addition  of  the  cross  sections  of  the  individual  elements.  It  can  also  be  seen  that  a  larger 
the  separation  between  the  ellipse  and  the  surface  results  in  an  increasingly  better 
approximation  by  the  incoherent  addition  with  respect  to  the  exact  result.  This  fact  has 
been  verified  by  examining  the  cumulative  root  mean  square  error  (cumulative  with 
respect  to  the  observation  angles). 


Figure  28:  6  wavelength  ellipse,  20  wavelengths  above  a  Gaussian  rough  surface 
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Figure  29:  6  wavelength  ellipse,  5  wavelengths  above  a  Gaussian  rough  surface 

Since  the  agreement  between  the  incoherent  addition  of  the  parts  and  the  exact  total 
cross  section  becomes  closer,  we  can  expect  a  threshold  for  the  distance  at  which  higher 
order  interactions  become  significant.  In  addition  to  the  separation,  this  threshold  will 
most  likely  depend  on  the  observation  direction,  the  illumination  direction  and  the 
orientation  of  the  ellipse.  Further  numerical  studies  will  be  required  to  find  these 
relationships.  We  can  see  that  the  calculation  of  the  returned  power  for  most  foliage 
components  will  probably  not  require  accounting  the  higher  order  interactions.  However, 
a  land-based  target  buried  beneath  the  foliage,  may  require  accounting  for  these 
interactions. 

4.3.5  Conclusions  and  Future  Efforts 

The  solution  of  the  cylinder  above  a  rough  surface  serves  as  a  basis  for  comparison 

with  the  first  order  multiple  scattering  approach  and  ultimately  the  radiative  transfer 
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approach.  The  results  from  these  simulations  will  justify  the  requirements  for  higher 
order  interactions  in  the  foliage-surface  scattering  problem.  As  one  would  expect,  we 
have  seen  an  increased  importance  for  the  higher  order  interactions  as  the  scatterer  is 
moved  closer  to  the  surface  or  its  size  increases.  This  method  is  obviously  inadequate  as 
the  number  of  scatterers  increases  or  as  the  incidence  angle  increases  since  the  problem 
becomes  numerically  intractable.  The  conclusions  drawn  from  the  example  are  not 
surprising  and  include  the  following 

•  the  iterations  required  increased  as  the  ratio  of  the  body  size  to  height  above  the 
surface  decreased 

•  the  interaction  between  the  surface  and  the  cylinder  decreases  as  the  separation  is 
increased  and  incidence  angle  is  decreased  (measured  from  the  vertical) 

•  the  beamwidth  must  be  significantly  larger  than  the  ellipse  to  see  the  effects  of 
multiple  scattering 

•  the  single/double  scatter  corrections  significantly  improve  the  estimate  for  power 
returned  relative  to  the  simple  power  addition  (at  least  for  small  roughness) 

In  addition,  we  have  seen  slower  convergence  in  our  technique  as  the  scatterer  is  moved 
closer  to  the  surface.  Consequently,  a  stabilized  bi-conjugate  gradient  solution 
(BiCSTAB)  in  combination  with  the  MOMI  has  been  implemented  as  an  aide  for  the 
convergence  of  the  problem.  We  note  that  the  application  of  the  BiCSTAB  routine  to  the 
MOM  equations  in  one  particular  example  (l  6k  cylinder  major  axis,  3k  above  the  rough 
surface)  did  not  converge  within  the  number  of  iterations  allotted.  Likewise,  the  straight 
MOMI  solution  required  30  iterations.  However,  when  the  BiCSTAB  routine  was  applied 
to  the  MOM  equations  after  the  MOMI  preconditioner  was  applied,  the  ellipse  and  rough 
surface  system  required  only  10  MOMI/BiCSTAB  iterations. 

Two  final  notes:  since  these  simulations  occur  with  monochromatic  waves  and  the 
interest  in  this  work  involves  pulsed  energy,  the  pulse  chosen  for  comparison  will  be 
slowly  varying  and  of  long  duration.  Primarily,  we  are  interested  in  the  importance  of  the 
interactions,  not  the  solution;  consequently,  until  a  time-dependent  code  is  introduced,  we 
will  assume  that  the  importance  of  the  interactions  in  the  pulsed  energy  problem  is 
similar  to  that  in  the  monochromatic  problem.  In  addition,  one  further  assumption  of  the 
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impulse  response  model  requires  our  attention;  the  assumption  of  no  wide-angle 
scattering.  This  may  be  significant  for  all  components  of  foliage  and  ground  based 
targets. 
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4.4  Approximate  Analytical  Solution  for  the  Moments  of  a  Single 
Scatterer  above  a  Rough  Surface 

If  it  is  found  that  the  scatterer  above  the  rough  surface  includes  important 
interactions  that  are  not  included  in  the  modified  first  order  multiple  scattering  solution, 
this  solution  must  be  refined.  In  addition,  if  the  return  from  a  strong  scattering  object, 
such  as  a  vehicle,  under  the  vegetation  is  desired,  then  the  coherent  return  may  be 
desired.  The  exact  solution  for  the  problem  of  multiple  scatterers  above  a  rough  surface 
will  become  numerically  intractable  as  the  number  of  scatterers  increase.  Consequently, 
alternative  methods  must  be  used  or  the  exact  solution  must  be  simplified.  As  we  have 
seen  in  the  previous  sections,  the  exact  solution  for  a  single  scatterer  above  the  rough 
surface  in  combination  with  the  first  order  multiple  scattering  will  produce  some  insight 
into  the  validity  of  the  radiative  transfer  result.  An  alternate  approach  that  simplifies  the 
exact  results,  yet  unlike  the  first  order  multiple  scattering  result,  maintains  the  coherent 
response,  begins  with  the  exact  integral  equations  and  incorporates  some  reasonable 
assumptions.  Like  the  first  order  multiple  scattering  response,  we  begin  with  a  single 
scatterer  over  a  rough  surface  and  propose  an  extension  to  N  scatterers  above  a  rough 
surface. 

We  start  with  coupled  integral  equations:  one  representing  the  current  on  the  rough 
surface  and  the  other  representing  that  on  the  scatterer.  From  equivalence,  the  MFIE  for 
the  current  on  the  scatterer  in  the  presence  of  the  rough  surface  can  be  written 

J„(fl)=2n1xHl(f,)+2n,x  JJJsl(r,')x  VG(r,^')dS, 

Scatterer 

Surface  f  ^ 

+  2n,x  JJJ,2(r2')*V'G(r„V)dS2  '  ' 

Rough 

Surface 

where  the  subscripts  1  and  '  2  will  indicate  points  or  currents  on  the  scatterer  and  the 
rough  surface,  respectively.  The  current  is  evaluated  at  the  observation  point,  which  is  on 
the  scatterer.  Note  the  presence  of  the  term,  the  last  term.  From  equivalence,  the  MFIE 
for  the  current  on  the  rough  surface  in  the  presence  of  the  scatterer  can  be  written 
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Ju(f2)=2n2  xfl'ft)+2n2  x  JJ  J„(f,')x  V'G(ri,f,')dS, 


Scatterer 

Surface 


+  2n2x  /JJs2(f2,)xV,G(r2,r^,)dS2 


(2) 


Rough 

Surface 


The  current  is  evaluated  at  the  observation  point,  which  is  on  the  scatterer.  Note  the 
presence  of  the  coupling,  the  second  term.  The  geometric  quantities  in  these  two 
equations  are  defined  in  Figure  30. 


Figure  30:  Geometry  for  the  Reduced  Integral  Equation  Approach 

4.4.1  The  Reduced  Integral  Representation 

The  overall  goal  is  to  simplify  the  solution  for  a  single  scatterer  over  a  rough 
surface.  Since  the  scatterer  is  assumed  to  be  small  with  respect  to  the  distance  to  the 
surface,  the  far-field  form  of  the  Green’s  function  will  be  used  for  interactions  involving 
the  scatterer  as  the  source  and  the  surface  as  the  observation  location.  In  addition,  when 
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the  object  is  modeled  as  a  smooth  ellipse  or  a  disc,  the  Physical  Optics  (PO) 
approximation  will  be  used  to  estimate  the  induced  currents  on  the  scatterer  due  to  the 
incident  field.  In  all  remaining  sections,  these  two  simplifications  shall  be  considered 
accurate.  Starting  with  the  coupled  integral  equations  and  substituting  into  the  integrands 
of  (1)  and  (2)  in  the  previous  section,  we  find  for  the  scatterer 

=  2n,xHl(r,)+2n1x  V'G(r„V)dS,  (1) 

Rough 

Surface 


In  this  formulation,  the  integral  equation  for  the  current  on  the  scatterer  still  involves  two 
unknowns.  For  the  currents  on  the  rough  surface,  we  find  the  following  equation 


Js2(r2)  =  2n2xH'(r2) 


xV'GffCr^r/jdS, 


+  2n2x  Jj Js2(r2')x V'G(f2,r2')dS2 

Rough 

Surface 


+  2n2  x 


II 


Scatterer 

Surface 


2n,  xH'(r,)  +  2n,  x  JJ  j„(V)x  V'G«^’)dS, 


Rough 

Surface 


(2) 

where  Gff(r,r’)  is  the  far-field  form  of  the  free  space  Greens  Function.  Those  vales  with 
a  subscripted  “1”  are  in  reference  to  the  scatterer  and  those  with  a  “2”  are  with  respect  to 
the  surface.  In  addition,  the  primed  coordinates  reference  the  sources  and  unprimed 
reference  the  observation  points.  This  integral  equation  consists  of  three  terms 

1 .  The  first  term  is  the  well-known  Kirchhoff  term 

2.  The  second  term  couples  the  currents  of  the  surface  to  that  of  the  scatterer 

3.  The  third  term  is  the  familiar  multiple  scattering  term  for  the  surface  to  surface 
interactions 

Notice  that  the  PO  current  on  the  scatterer  is  known;  consequently,  the  only  unknown  in 
the  integral  equation  for  the  surface  current  is  the  surface  current  itself.  Moving  the 
normal  unit  vectors  inside  the  integrals  and  using  the  vector  identity 

A  x  B  x  C  =  B(AC)-C(A-B) 
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(3) 


Expanding  the  gradient  of  the  Greens  Function: 


V'G(r.r')  =  ^r4r  +  jk  G(r =  — L-  +  jk  G(r,f)k,  (4) 
r-r'  r-r  r-  r' 


where,  following  the  notation  of  [Ishimaru,  1978],  the  direction  from  the  source  point  to 
the  observation  point  is  given  by 

a  (r  ~  r')  .  . 

k  =  -j - -  =  the  scattering  direction 

|r  -  r’J 

Employing  the  far-field  approximation,  the  Greens  function  and  its  gradient  become 


G(r,r’)  =  ~r"e’jfel  s  Gff(r,r') 
47tR 


V'  G(r,r')  =  jki— e-j^ks  =jkG(r,r')ks 
4nR 


where  R  s  r-r'.  Further  reduction  of  integral  equations  is  accomplished  by  assuming 
the  scatterer  has  a  definite  geometrical  shape  (disc,  etc.).  In  this  case  the  backscatter  and 
forward  scatter  from  known  scatterer  geometry  may  have  an  analytical  result,  further 
simplifying  the  integral  equations. 


4.4.2  Reduction  for  a  Circular  Disk  (3-D)  Scatterer  above  a  Rough  Surface 

In  3-D,  the  integral  equations  were  specialized  to  a  circular  disc  with  random  tilt  and 
height  above  a  randomly  rough  surface  and  in  2-D,  the  integral  equations  were 
specialized  to  a  strip  with  random  tilt  above  a  corrugated  surface.  Only  the  results  of 
these  derivations  will  be  given  in  this  report.  These  results  should  include  more 
interaction  terms  with  the  surface  than  the  first  order  multiple  scattering  theory  but  with 
less  computational  demand  than  the  exact  solution.  For  a  flat  disc,  horizontally  suspended 
over  a  rough  surface,  the  current  on  the  rough  surface  is 
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Js(r2)  =2n2xHi(f2)  +  jk  CXpj  j^25  ^  exp{-  j  kR2S  •  zh} 

2ttR2S 

'  //[(“  ’  R2S )(z x  H s  (r, ' ' ) )  +  R2S  (n2  •  [z  x  fi5  (r, ' )  ]  )]exp{-  j  kR2S  •  r,  '}dS, ' 

si 

+  2  JJ{ Js (%' ■ )  (n,  •  RJr )  - (n2  •  Js (f2 '  ))ft  „  }f 9' G(f; ,  r2 ' ' )| dS2 1 ' 

S2 


+  eXphjkR-s}kAexp 


4tc  R 


2s 


f-jkh: 

1  R2S 


Ji(kA) 


•R2S)  (^2  ‘  Js(r2  ))R2s]-d  [(n2  •  R2S)rs2  (^2 

Xvc 


^S2 


exp{- jkRS2} _ f-jkfa 


4tt2R 


exp< 


s2 


R 


•dS, 


s2 


(1) 


where  the  Bessel  Function  of  the  first  kind,  Ji(x),  has  arisen  due  to  the  circular  disk  and  A 
is  the  area  of  one  face  of  the  disk.  It  has  been  assumed  that  the  disk  is  very  thin.  The  first 
term  is  the  Kirchhoff  current  on  the  surface.  The  second  term,  the  first  integral  term,  is 
an  additional  Kirchhoff  current  term  due  to  the  incident  field  from  the  Kirchhoff  current 
on  the  disc:  the  Kirchhoff  current  on  the  disk  radiates  to  the  surface.  Shadowing  must  be 
accounted  for  in  the  use  of  this  term:  the  Physical  optics  current  on  the  underside  of  the 
disc,  due  to  the  incident  field  will  typically  be  zero.  The  third  term  is  the  surface  in 
isolation.  Finally,  the  fourth  term  is  a  multiple  interaction  result:  the  current  on  the 
surface  radiates  to  the  disk  and  is  re-radiated  to  the  surface  again  ...  ad  nauseum.  Note 
that  this  integral  equation  has  only  one  unknown:  the  current  on  the  surface. 


The  integral  equation  for  the  disk  simply  involves  Kirchhoff  current  and  the  current 
due  to  the  surface  radiating  to  the  disk.  Note  that  this  equation  is  fully  coupled  to  the 
solution  for  the  current  on  the  surface 
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If  the  surface  is  gently  undulating,  the  z-directed  currents  will  be  nearly  zero,  hence,  the 
integral  equations  for  the  currents  become 

Js(^)  =2n2  x Hj(r2)  +  jk  £X~j— pk~S'"exp{-  jkR2S  •  zh} 


'  jj[(-n2  •R2S)(zxHj(fi’))+  R2S(h2-[zxHi(f1,)])]exp{-jkR2S-fI’}dS1' 

SI 

+  2  JJ{ Js(r2')  (h2  •R22.)-(n2-Js(f2,))R22.  }|V'G(f2,r2')|dS2’ 

S2 


+ 


exp{-  jkR2S} 
4Tt2R2S 


J,(kA) 


J,(f,)=2zx  H,(?,)  +  2JJJS!(V) 

S2 


1  }  h 

+  ik  3 -rrrGfe, v)  ^ 

r."r2  J  r!  -  r2  | 


(4) 

A  more  complex  result  is  available  for  a  disk  with  arbitrary  tilt  and  in  2-D,  a  strip  with 
arbitrary  tilt. 
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4.4.3  Conclusions  and  Future  Work 

This  section  has  started  with  the  exact  coupled  integral  equations  for  a  rough  surface 
in  the  presence  of  a  single  scatterer  and  reduced  these  using  the  far-field  form  of  the 
Green’s  function  for  the  scattered  field  due  to  the  scatterer  and  the  physical  optics 
solution  to  the  scatterer’ s  current.  Furthermore,  we  have  isolated  an  approximate  integral 
equation  for  the  surface  scattering  in  the  presence  of  the  disk.  This  integral  equation  only 
involves  one  unknown,  the  current  on  the  rough  surface. 

In  addition  we  have  simplified  the  expression  for  a  scatterer  that  is  a  circular  disk.  In 
addition,  the  term  that  represents  the  multiple  interactions  between  the  scatterer  and  the 
surface  has  been  isolated  and  should  be  evaluated  relative  to  the  surface  in  isolation.  This 
investigation  will  yield  an  analytical  solution  as  to  the  validity  of  the  single  interaction 
assumption  for  the  radiative  transfer  result.  The  next  step  in  this  process  is  to  numerically 
implement  the  above  integral  equations.  These  results  may  then  be  compared  with  those 
obtained  with  MOM/MOMI.  Concurrently,  the  average  results  will  be  attempted 
analytically  for  a  random  height  and  then  the  random  orientation.  These  results  should  be 
implemented  numerically. 

In  the  extension  to  N  scatterers  above  a  rough  surface,  no  interaction  between 
scatterers  will  be  accommodated;  the  result  will  include  interaction  with  the  surface,  like 
the  first  order  multiple  scattering  result.  Unlike  the  first  order  multiple  scattering  result 
discussed  earlier,  these  results  will  include  a  more  comprehensive  treatment  of  the 
interaction  with  the  surface  in  addition  to  the  preservation  of  the  coherent  field. 
Numerical  solutions  and  comparison  with  MOM  results  for  two  to  three  scatterers  can  be 
performed  in  an  effort  to  assess  the  mutual  interactions  among  the  scatterers  themselves 
(ignored  by  the  presented  single  scatter  theory).  An  attempt  to  derive  analytic  expressions 
for  the  field  moments  from  these  equations  will  be  made,  including  the  mutual  coherence 
functions. 

In  addition  to  a  volume  return  component,  this  reduced  integral  equation  formulation 
result  will  produce  the  most  comprehensive  treatment  of  the  scattered  field  from  a 
collection  of  scatterers  above  a  rough  surface.  We  note  that  many  useful  methods  already 
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exist  which  predict  both  the  coherent  and  incoherent  responses  from  a  volume  of  discrete 
scatterers  that  may  be  combined  with  these  results  in  order  to  produce  a  comprehensive 
model.  These  include  the  DWBA  [Lang,  1981],  Cumulative  Forward  Scatter,  Single 
Backscatter  theory  [deWolf  ,1971],  its  extension  to  pulse  propagation  [Ishimaru,1980] 
and  the  related  spectral  approach  [Rino,  1988]. 
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5  Conclusions  and  Future  Action 


In  predicting  the  radar  return  from  vegetation,  a  number  of  approaches  have  been 
developed.  Typically,  the  overly  cumbersome  radiative  transfer  seems  to  be  very  popular. 
Exceptions  to  this  trend  have  been  found  in  the  works  of  [Schwering,  1985],  and  [Brown 
and  Adams,  1998a],  Schwering  has  used  a  simple  scattering  function  to  represent  the 
scattering  properties  of  all  components  of  the  volume.  A  notable  wave  approach  to  this 
problem  is  that  of  Lang  [Lang,  1981]  who  has  used  the  Distorted  Wave  Bom 
Approximation.  In  this  study  we  have  formulated  a  simple  model  that  is  numerically 
efficient  and  depends  on  the  empirical  identification  of  effective  parameters.  However, 
there  are  a  number  of  verification  studies  that  must  be  performed  and  several  different 
levels  of  verification  have  been  outlined  in  this  report.  The  major  thrust  of  these 
verifications  is  to  identify  the  necessary  level  of  approximation  for  the  foliage-surface 
interaction. 

The  Impulse  Response  (Convolutional)  model  allows  a  superposition  of  surface  and 
volume  responses.  Its  numerical  implementation  is  via  the  Fast  Fourier  transform,  FFT, 
which  allows  a  fast  numerical  solution.  This  approach,  although  originally  derived  from 
the  radar  equation  and  then  by  the  radiative  transfer  equations,  has  been  shown  to  be 
equivalent  to  the  “first  order  multiple  scattering  theory.”  This  equivalence  has  been 
derived  under  the  assumptions  of  narrow-band  and  narrow-beamwidth  with  a  limited 
scattering  pattern.  From  this  equivalence,  we  have  proposed  a  method  using  the  first 
order  multiple  scattering  to  verify  the  single  passage  assumption  for  foliage-surface 
interaction  that  is  inherent  in  the  radiative  transfer  approach. 

The  future  direction  of  this  work  will  involve  thoroughly  investigating  the 
assumptions  of  the  impulse  responses  and  consequently,  establishing  the  range  of  their 
validity;  these  investigations  may  in  turn  expand  the  applicability  of  the  model. 
Generalizing  the  model  to  include  lower  frequencies,  larger  beamwidths,  etc.  will  require 
an  investigation  into  the  range  of  validity  of  each  of  the  above  limitations  and  may 
involve  generalization  of  the  model.  Some  of  these  investigations  and  extensions  have 
been  performed.  The  most  straightforward  model  improvement  will  be  the  addition  of 
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polarization.  The  use  of  the  Stokes  vector  will  allow  for  polarization  dependence  in  the 
model.  Depolarization  will  enter  naturally  through  the  Stokes  vector  definition,  existing 
models,  and  measured  data.  After  accounting  for  the  polarization  effects,  the  incident 
beam  waveform  can  be  given  a  pulse  shape.  This  pulse  will  be  constructed  using  Fourier 
analysis  in  order  to  recreate  pulse  scattering  with  an  incident  beam  using  a  full  radiative 
transfer  approach.  This  pulsed  waveform  result  is  then  directly  comparable  to  our 
impulse  response  model.  For  large  bandwidth  signal,  each  Fourier  component  of  the 
pulse  may  be  treated  separately  in  radiative  transfer  theory;  hence,  the  scattering 
characteristic  of  the  scatterers  can  be  changed  with  frequency.  When  derived  from  the 
single  (or  higher  order  multiple)  scattering  approach,  the  impulse  response  model  will 
require  a  more  general  form  for  the  two  frequency  mutual  coherence  function. 
Investigations  into  the  other  approximations  will  require  investigations  into  the  effects  of 
wide-angle  and  multiple  scattering. 

The  effects  of  multiple  scattering  may  be  established  by  inserting  the  next  higher- 
order  scattering  correction  in  our  model.  Through  these  corrections,  we  may  estimate  the 
system  parameters  or  foliage  and  surface  conditions  that  will  either  cause  the  model  to 
fail  completely,  e.g.  become  grossly  invalid,  or  give  rise  to  significant  changes  in  the 
expected  effective  parameters  in  the  model.  Beginning  with  the  development  of  the 
impulse  response  method  using  multiple  scatter  theory,  we  have  developed  a  rigorous 
formalism  that  will  simplify  to  the  impulse  response  model  under  the  given  assumptions. 
Under  this  model,  higher  order  scattering  corrections  may  be  added  in  each  of  the 
interaction  regions:  foliage-foliage,  surface-surface  or  foliage-surface.  For  example, 
known  numerical  techniques  may  be  applied  to  correct  for  multiple  scattering  among 
surface  elements;  the  applicability  of  the  facet  model  used  by  the  impulse  response  model 
may  be  directly  compared  with  these  numerical  techniques.  The  multiple  interactions  and 
the  effect  of  wide-angle  scattering  between  the  foliage  constituents  and  the  surface  can 
also  be  implemented  numerically  but  this  creates  a  numerically  intractable  problem  as  the 
number  of  foliage  components  grows.  Hence,  we  have  looked  to  reduce  this  problem  and 
only  characterize  the  important  interactions  between  the  foliage  and  the  surface. 
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Although  the  mutual  interactions  between  the  scatters  comprising  the  vegetation  are 
neglected  in  our  first  order  model,  some  interaction  between  each  of  these  scatterers  and 
the  surface  may  be  non-negligible.  In  the  first  order  model,  the  backscattered  power 
density  from  the  rough  surface  and  the  foliage  volume  separate  into  two  distinct 
components.  Consequently,  by  construction,  this  model  assumes  that  there  is  no 
interaction  between  the  surface  and  the  volume  scattering  elements  other  than  the 
exponential  decay  of  the  surface  scattered  power  through  the  volume.  In  the  course  of  our 
research,  we  have  developed  an  exact  integral  equation  approach  that  simulates  the 
interaction  of  a  one-dimensional,  statistically  rough  surface  in  the  presence  of  an 
elliptical  cylinder  at  some  distance  above  the  surface.  This  effort  is  necessary  in  order  to 
compare  the  first  order  multiple  scattering  result  with  a  result  with  known  accuracy. 
These  results  provide  a  basis  of  comparison  for  the  impulse  response. 

Implemented  via  the  MOMI  method  of  solution  to  the  MOM  problem,  this  numerical 
model  permits  the  specification  of  each  order  of  interaction  that  is  present  in  the  result, 
e.g.  the  first  and  second  order  interactions,  between  the  scatterer  and  the  surface.  To 
accommodate  strong  interactions  between  surface  and  foliage,  several  advanced 
numerical  techniques  have  been  (will  be)  implemented,  e.g.  the  method  of  multiple 
ordered  interactions,  the  stabilized  bi-conjugate  gradient  algorithm  and  the  fast  multipole 
method.  The  results  presented  are  restricted  to  the  TE  incidence  but  they  have  also  been 
extended  to  the  construction  of  the  MOMI  model  for  TM  polarization  for  an  elliptical 
cylinder  over  a  rough  surface.  We  have  found  that  a  simple  incoherent  addition  of  the 
scattered  power  will  be  accurate  for  most  foliage  components;  this  will  depend  on  the 
ratio  of  the  scatterer’ s  largest  dimension  to  its  displacement  above  the  rough  surface.  A 
large  object  such  as  a  target  embedded  within  the  foliage,  on  the  other  hand,  will  require 
multiple  scattering  and  wide-angle  corrections  for  an  accurate  scattering  prediction.  It  has 
been  observed  that  the  interaction  between  the  surface  and  a  nearby  object  can  be 
significant,  particularly  when  the  scatterer  is  closer  to  the  surface  and  as  incidence  angle 
is  further  removed  from  nadir. 

If  these  interactions  between  surface  and  scatterer  have  a  significant  effect  on  the 
returned  power  waveform,  we  may  resort  to  the  full  first  order  multiple  scatter  theory. 
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which  includes  wide-angle  scattering.  This  theory  will  more  fully  account  for  foliage- 
surface  interaction  and  extend  the  impulse  response  model.  Although  this  will  add 
complexity  to  the  model,  some  advantages  of  the  convolutional  approach  may  be 
retained.  A  plan  has  been  presented  to  extend  the  first  order  multiple  scattering  solution 
for  foliage  and  foliage-surface  interaction  to  both  broadbeam  and  broadband  cases:  this  is 
not  possible  in  the  radiative  transfer  approach  as  presented.  By  including  double 
scattering  events  (scattering  between  the  foliage  and  every  surface  component),  the 
assumption  of  pure  forward  scatter  and  backscatter  will  be  abandoned.  Through  this 
extension  of  the  impulse  response  method  each  component  of  the  foliage  will  interact 
with  every  surface  element  (facet)  rather  than  a  single  surface  element.  This  numerical 
model  may  also  be  used  to  establish  the  multiple  scattering  significance  for  a  target 
embedded  within  our  vegetated  surface.  Such  a  target  may  be  found  through  a 
comparison  of  the  return  waveform  from  the  vegetated  rough  surface  with  that  of  the 
vegetated  rough  surface  with  an  embedded  target.  We  may  then  tailor  our  system 
parameters  (pulse  length,  beamwidth,  etc.)  to  more  readily  identify  its  presence. 

The  full,  first-order  multiple  scatter  theory  will  also  test  the  validity  of  the  narrow 
beamwidth  assumption  in  the  foliage-surface  interaction.  However,  for  the  foliage-foliage 
interaction  another  more  complete  test  will  involve  a  solution  of  the  unsimplified 
radiative  transfer  equation.  This  solution  will  test  the  effects  of  higher  orders  of  multiple 
interactions  on  the  received  waveform.  This  assumption  presumes  that  the  returned  power 
is  confined  to  a  small  region  around  the  transmitting  beam  axis.  Beam  broadening 
through  wide-angle  scattering  and  multiple  scattering  in  a  tenuous  media  can  be 
neglected  if  the  receiver  and  transmitter  beamwidths  are  narrow  since  multiple  scattering 
effects  and  wide-angle  scattering  are  negligible  within  a  small  illuminated  volume. 
Through  comparison  with  the  impulse  response  model,  a  range  of  validity  for  the  narrow 
beam  assumption  will  be  established.  Although  the  steady  state  solution  of  the  radiative 
transfer  equations  for  a  plane  wave  incident  to  a  random  media  is  widely  available  in 
literature,  this  solution  must  be  extended  to  account  for  an  actual  antenna  pattern.  To  this 
end,  we  have  constructed  a  radiative  transfer  solution  for  beam  incidence  and  are  in  the 
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process  of  evaluating  the  role  of  multiple  scattering,  wide-angle  scattering  and  beam 
broadening. 

In  addition  to  the  impulse  response  approach,  a  second  analytical  approach  has  been 
formulated  which  begins  with  the  exact  integral  equations.  We  hope  to  find  an  analytical, 
or  at  least  computationally  efficient,  solution  to  both  the  mean  and  second  order,  time 
dependent  power  density  for  the  single  scatterer  above  a  rough  surface.  Like  the  impulse 
response  approach,  the  incoherent  results  for  the  superposition  of  N  scatterers  will 
provide  a  new  look  at  the  single  scattering  theory  for  an  ensemble  of  particles  above  a 
rough  surface  and  should  serve  as  a  check  against  the  impulse  response  result.  However, 
unlike  the  impulse  response  method,  this  approach  will  provide  the  coherent  field 
returned  from  the  volume  and  the  rough  surface  and  it  will  include  higher  orders  of 
multiple  scattering.  This  result  should  be  directly  comparable  to  the  Distorted  Wave  Bom 
Result  (DWBA),  see  [Lang,  1981].  Hence,  one  additional  task  includes  constructing  a 
DWBA  model  including  a  rough  surface,  based  on  Lang’s  work  and  in  order  to  compare 
with  this  reduced  integral  equation  result. 
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Appendix  A  First  Order  Multiple  Scattering 

The  simplest  approach  to  propagation  through  a  random  media  is  via  single  scatter 
theory.  Single  scatter  theory  retains  its  simplicity  by  incorporating  the  assumption  that 
the  field  incident  on  a  random  media  interacts  with  each  scatterer  only  once  and  no 
multiple  interactions  occur  among  the  scatterers  that  encompass  the  random  media.  This 
creates  a  scattered  field,  which  is  a  power-like  summation  of  scattered  power  from  each 
scatterer. 

In  a  previous  section,  the  impulse  response  approach  was  introduced  as  an  efficient 
means  for  computing  the  average  intensity  from  a  rough  surface.  In  this  section  this 
impulse  response  technique  is  again  extended  to  a  tenuous  (sparse)  random  media 
covering  a  rough  surface.  Like  the  radiative  transfer  extension,  this  method  will  use  the 
convolutional  approach.  Unlike  the  extension  from  the  radiative  transfer  theory,  the 
following  method  applies  the  convolutional  approach  to  a  strongly  scattering,  yet  tenuous 
media,  which  accounts  for  scattering  out  of  the  radial  path  with  respect  to  the  antenna;  the 
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radiative  transfer  approach  only  accounts  for  strictly  forward  scattering  and 
backscattering. 


A.l  Review  of  Classical  Single  Scatter  Theory 

Assuming  plane  wave  is  incident  on  the  scatterer,  the  incident  field  and  the  associated 
mean  power  density  in  free  space  are  written  as 

E,(r)=  Eoe-^-'ki,  S,  =  i  (§,(?)  x  H*(?))= 


(1) 

A 

where  the  incident  direction  is  denoted  by  ks  and  the  impedance  of  free  space  is  given  as 

“Ho  = 

In  the  development  of  the  single  scatter  theory  in  a  random  media,  the  total  received 
power,  PR,  is  written  as  a  summation  of  the  power  scattered  once  by  each  particle.  This 
simplification  is  due  to  the  randomness  of  the  scatterers  in  the  media:  all  interference 
effects  are  neglected  [Ishimaru,  1997].  The  scattered  field,  Es(r),  and  power  density,  Ss, 
due  to  the  scatterer  can  be  represented  as 


Es(f)  =  ?(£„£,  )eXp(RjkR)  .  Ss  =i(Es(f)xH;(r))=^-|Es(rfks 


where  f(ks,k;)  is  the  scattering  amplitude,  R  is  the  range  from  the  radar  to  the  particle 


A  A  • 

and  the  incident  direction  is  denoted  by  k;  and  the  scattering  direction  is  denoted  by  ks . 
Next  we  make  several  definitions  concerning  the  radar  cross-section  of  this  scatterer. 
First  is  the  differential  radar  cross-section,  ad  for  given  incident  and  scattered  field 
directions 
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(3) 


Consequently,  the  backscattering  cross-section  will  be  denoted  as  ab  which  will  imply 
the  following:  47iad(ki,-ki).  Another  important  cross-section  in  this  work  is  the  total 
observed  cross-section  (the  scattering  cross-section),  cs .  Denoting  a  differential  unit  of 
solid  angle  as  dQ,  this  cross-section  is  written 


cjs=  ffad  dQ 
4n 


In  the  development  of  the  single  scatter  theory  in  a  random  media,  the  total  received 
power,  PR,  is  written  as  a  summation  of  the  power  scattered  once  by  each  particle.  This 
simplification  is  due  to  the  randomness  of  the  scatterers  in  the  media:  all  interference 
effects  are  neglected  [Ishimaru,  1997].  Consequently,  summing  up  the  power,  PR, 
returned  from  a  continuum  of  scatterers  to  the  radar  from  a  random  media  due  to  the 
transmitted  power,  PT,  can  be  expressed  using  the  standard  radar  equation 


PR  =P 


.  x2  [G(e,  4.)] 2  P  crb  (e,  <|>)  dy 

V 


(4tt)3R4 


where:  X  =  wavelength  of  the  carrier 

G(0,  <j>)  =  radar  antenna  gain  in  the  direction  (0,<j>)  or  kj 

<7b(0,  <j>)  =  particle  backscattering  cross  section  per  unit  area 
dV  =  elemental  volume 

R  =  slant  range  from  the  radar  to  dV 
p  =  particle  density  per  unit  volume 


(4) 


In  this  expression  the  transmitted  power  is  assumed  to  be  a  continuous  wave  signal; 
consequently,  the  frequency  dependence  is  monochromatic  and  monochromatic 
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dependence  is  implicit  in  the  backscattering  cross-section  and  the  transmitting  and 
receiving  subsystems  as  well  as  any  assumptions  related  to  the  intervening  media. 


In  an  attempt  to  improve  this  formulation,  an  accounting  for  the  loss  into  the  media 
from  the  source  to  the  scatterer  has  results  in  first  order  multiple  scattering  and  takes  the 
a  modified  form  with  respect  to  the  previous  expression  for  the  received  power  [Ishimaru, 
1997].  This  modified  form  can  be  described  as  follows.  Once  the  absorption  of  the 
scatterer  is  included,  this  lost  energy,  represented  by  the  absorption  cross  section,  aa ,  is 
added  to  the  scattering  cross  section  to  form  the  total  or  extinction  cross  section, 
ot  =  ct3  +  cts.  Hence,  as  a  coherent  wave  travels  through  an  uncorrelated  random 
media,  the  coherent  field  is  diminished  by  the  total  cross  section  of  the  scatterers 
encountered.  Hence  over  a  differential  change  in  distance,  the  coherent  power  density  can 
be  written  as  in  the  following  expression 


Hence  the  incident  power  density  at  the  elemental  scattering  area,  dV,  at  a  depth  R  into 
the  media  can  be  written  as 

Sj (R; kj )  =  Si(R0;ki)exp{-  J^p(at)dR} 


where  (a, )  is  the  average  total  cross  section.  In  terms  of  the  transmitted  power 

PT  G(0, <J>)  exp{-  £  po,dR  } 

p.  _ 

4tiRz 

Note  the  lack  of  any  forcing  or  source  terms  since  energy  cannot  scatter  into  the  coherent 
field  when  the  media  is  random  with  no  correlation  among  the  particles.  The  power 
received  by  the  antenna  is  then  expressed  as 
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Pr  -  PT 


JH 

volume  of 
scatterers 


p{-2J0V,dR}dV 

(47r)3R4 


(5) 

This  power  density  returned  is  identical  to  the  power  density  predicted  by  the 
convolutional  radiative  transfer  result;  however,  the  integrations  have  been  re-arranged 
and  time  dependence  is  ignored  at  this  point. 
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A.2  The  Scattered  Pulse 


In  order  to  construct  the  response  of  the  media  to  a  scattered  pulse,  the  correlation  of 
the  output-scattered  fields  must  be  derived.  Loosely  following  the  notation  and 
development  of  Ishimaru  [1997],  a  general  expression  for  the  correlation  is  found  in  the 
following  equation 

Bu(t,,t2)  =  fdaj  dcD2Uj(co,)U*(co2)r exp{- j^t,  -co2t2)} 

•*—00  CO  1 


(6) 

where  U;  (co)  is  the  complex  envelope  of  the  incident  wave  form  at  the  time  harmonic 
frequency  co  and  F  is  the  two-frequency  mutual  coherence  function.  The  two-frequency 
mutual  coherence  function  is  the  correlation  of  the  time-varying,  frequency  domain 
transfer  function,  H((o,t),  at  two  different  frequencies  and  two  different  times 

r  =  r(co0  +  +<o2;t„t2)  =  (H(co0  +  +  a2,t2)) 

(7) 

Once  the  two-frequency  mutual  coherence  is  constructed,  the  scattered  intensity  is  found 
when  ti  =  t2  =  t  and  r  ->  rQ  [Ishimaru,  1997] 

I(t)=f  dcoif  d<02Ui(©i)Ui’(©2)ro  exp{- j(©i  -©2)t} 


(8) 

In  first  order  multiple  scattering  theory,  the  transfer  function  at  a  single  frequency  can  be 
simply  derived  from  the  transform  of  the  spatial  scattering  function  f(ks,k;)  modified  by 
the  loss  in  the  media  and  gain  of  the  antenna  at  the  given  frequency.  For  simplicity, 
antenna  gain  at  a  given  frequency  for  monostatic  operation  can  be  represented  as  follows 

GR(e,<J>;co)  =  GT(0, <{>;©)  =  G(co) 


(9) 
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Note  the  monostatic  assumption:  transmitting  and  receiving  directions  are  identical. 
Next,  we  rewrite  the  gain  of  the  transmit  antenna  as 


_1_ 

4n 


GT(ki,fi>)  =  gT(ki5co) 


(10) 


Additionally,  assume  that  the  gain  on  transmit  and  the  effective  aperture  area  on  receive 
are  represented  by 


"7“  GR  (oj)  =  g*  (to) 
471 


(11) 


Assuming  the  media  is  stationary  in  time,  the  frequency  domain  transfer  function  at  the 
single  frequency  co,  can  be  expressed  as  follows  [Ishimaru,  1 997] 


-2jkR 


H(o)  =  gT  (o)gR  (co)  F(co)  exp(-  2  £  p(ot)dR 


(12) 


where,  the  Fourier  transform  of  the  scattering  function  at  a  given  frequency  yields 

F(g>)  =  F(©;ks,kj)  =  3{f(kf,ks) } 

Once  the  two-frequency  mutual  coherence  function  is  constructed,  the  backscattered 
intensity  has  been  found 


volume  of 
scatterers 


I(t)=  IJJ  |  jA^.t-^Je^dco,  Ui.,1-  c  } 


2R 


eJO)2tdco, 


pdV 


(13) 


where 


-2jkR 


A(t,on)  =  Ui(con)F(con;ks,ki)gT(co)gR(co)-^1-expj-2|)  p(at)dR 
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Uj(con)  =  the  complex  envelope  of  the  incident  signal  at  the  frequency,  con 

and  the  asterisk  signifies  the  complex  conjugate.  When  the  bandwidth  of  the  pulse  is 
narrow  with  respect  to  the  carrier  frequency,  the  narrow-band  approximation  can  be 
made.  In  this  case,  the  scattering  function,  as  a  function  of  frequency,  is  roughly  constant 
and  can  be  evaluated  at  the  center,  carrier  frequency.  Once  this  narrow  bandwidth 
approximation  assumption  is  made,  the  only  frequency  dependence  in  the  two-frequency 
mutual  coherence  function  appears  in  the  Fourier  kernel.  A  change  of  variables  to  the 
difference  frequency  ©d  =  <»[  -  co2  yields  a  simpler  expression.  The  two  frequency  mutual 
function  for  a  narrow  band  input  signal  becomes  [Ishimaru,  1 997] 

r.  HFK)!’ 

K 

(14) 

When  a  finite  number  of  discrete  scatterers  is  present,  in  contrast  to  the  formulation 
above,  the  volume  integration  over  the  density  of  scatterers  will  be  replaced  with  a 
discrete  summation. 

With  the  narrow  beamwidth  and  narrow  band  approximations,  the  inverse  transforms 
are  easily  performed  and  the  backscattered  waveform  reduces  to  [Ishimaru,  1997] 

1(0  =  JIJ  ft  P  gb  (9>  40  exp { -  2 J  *  p(o , )dR  }  u,f 1  -  — )  dV 

volume  of  (471 )  R  \  C  J 

scatterers 

=  SIS  f 1  p  <t„  (9. 40  exp{-  21 1  k.  (R)dR  }  PT  ft  -  —  ]  dV 

volume  of  (471 )  R  ^  C  ) 

scatterers 


(15) 

Rearranging  the  order  of  integration,  this  can  be  seen  to  be  equivalent  to  the  radiative 
transfer  result  derived  in  Chapter  3.  Assuming  a  the  speed  of  light  is  the  same 
everywhere,  the  radiative  transfer  approach  yields 


I(-r;r,9,<j>,t) 


f  r,0+5(x)sec6+d,  sec9 
Jr10+£(x)secei 


cb(a)I 


104 


When  we  integrate  this  radiative  transfer  result  over  the  surface,  substitute  for  the 
incident  intensity,  simplify  the  limits  and  compensate  for  the  effective  receiving  area  of 
the  antenna,  identical  results  are  realized  for  the  radiative  transfer  (see  Chapter  3)  and  the 
first  order  multiple  scattering  approaches. 


oo  2x 


^  G  (0,jO  p 

4  rT 


t  - 


2a 

c 


o  j 


exp  |  -  2  ke  (p)  dp  j  da  pd(j)dp 


(16) 

In  addition,  note  that  the  particle  density  per  unit  volume,  p,  is  assumed  to  be  unity. 

Like  the  radiative  transfer  result,  the  first  order  multiple  scattering  result  can 
accommodate  a  spatially  varying  velocity  and  may  be  re-cast  into  a  convolutional, 
impulse  response  form  when  the  narrow-band  and  narrow-beamwidth  approximations  are 
employed.  More  importantly,  a  further  limitation  of  the  radiative  transfer  approach  has 
been  identified:  use  of  a  narrow  bandwidth  approximation.  This  assumption  is  expected 
due  to  the  use  of  constant  forward  and  backscatter  coefficients  with  respect  to  frequency 
in  deriving  the  radiative  transfer  results.  The  use  of  the  narrow  beamwidth  approximation 
has  already  been  identified  in  the  radiative  transfer  approach.  However,  using  the  full 
expression  for  two-frequency  mutual  coherence  function,  it  is  possible  that  the  impulse 
response  approach  may  be  extended  to  broader  bandwidth  pulses  in  addition  to  broader 
beamwidth  antenna  patterns  while  maintaining  some  convolutional  aspects.  This  premise 
is  still  under  investigation. 
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ABSTRACT 

We  derive  the  nonhomogeneity  detector  (NHD)  for  non- 
Gaussian  interference  scenarios  and  present  a  statistical  anal¬ 
ysis  of  the  method.  The  non-Gaussian  interference  scenario 
is  assumed  to  be  modeled  by  a  spherically  invariant  random 
process  (SIRP).  We  present  two  methods  for  selecting  rep¬ 
resentative  (homogeneous)  training  data  based  on  our  sta¬ 
tistical  analysis  of  the  NHD  for  finite  sample  support  used 
in  covariance  estimation.  In  particular,  exact  theoretical  ex¬ 
pressions  for  the  NHD  test  statistic  probability  density  func¬ 
tion  (PDF)  and  its  moments  are  derived.  Additionally,  we 
note  that  for  SIRP  interference,  a  simple  transformation  of 
the  NHD  test  statistic  admits  an  elegant  representation  as 
the  ratio  of  a  central-F  distributed  random  variable  and  a 
beta  distributed  loss  factor  random  variable. 

1.  INTRODUCTION 

An  important  issue  in  space-time  adaptive  processing  (STAP) 
for  radar  target  detection  is  the  formation  and  inversion  of 
the  covariance  matrix  underlying  the  clutter/  interference. 
In  practice,  the  unknown  interference  covariance  matrix  is 
estimated  from  a  set  of  independent  identically  distributed 
(iid)  target-free  training  data  which  is  assumed  to  be  rep¬ 
resentative  of  the  interference  statistics  in  a  cell  under  test. 
Frequently,  the  training  data  is  subject  to  contamination  by 
discrete  scatterers  or  interfering  targets.  In  either  event, 
the  training  data  becomes  nonhomogeneous.  As  a  result, 
it  is  non  representative  of  the  interference  in  the  test  cell. 
Estimates  of  the  covariance  matrix  from  nonhomogeneous 
training  data  result  in  severely  undemulled  clutter.  Conse¬ 
quently,  CFAR  and  detection  performance  suffer.  Signifi¬ 
cant  performance  improvement  can  be  achieved  by  employ¬ 
ing  pre-processing  to  select  representative  training  data. 

The  problem  of  target  detection  using  improved  training 


strategies  has  been  considered  in  [1-8],  The  impact  of  non¬ 
homogeneity  on  STAP  performance  is  considered  in  [8-1 1]. 
The  works  of  [1-4, 8, 12]  have  addressed  the  use  of  the  non¬ 
homogeneity  detector  (NHD)  based  on  the  generalized  in¬ 
ner  product  (GIP)  measure  for  STAP  problems  involving 
Gaussian  interference  scenarios.  This  work  was  extended 
significantly  in  [13, 14]  to  include  the  effects  of  finite  sam¬ 
ple  support  used  for  covariance  matrix  estimation.  How¬ 
ever,  the  corresponding  problem  for  non-Gaussian  interfer¬ 
ence  scenarios  has  received  limited  attention. 

In  this  paper,  we  derive  the  NHD  for  non-Gaussian  in¬ 
terference  scenarios,  which  can  be  modeled  by  spherically 
invariant  random  processes  (SIRP)  and  present  a  statistical 
analysis  of  the  resultant  NHD  test  statistic.  In  general,  the 
problem  of  non-homogeneity  detection  for  SIRPs  is  quite 
difficult  due  to  the  fact  that  the  underlying  SIRP  covari¬ 
ance  matrix  and  characteristic  PDF  are  unknown.  For  con¬ 
venience,  knowledge  of  the  SIRP  characteristic  PDF  is  as¬ 
sumed  in  this  paper. 

2.  PRELIMINARIES 

Let  x  =  [aq  x2  ...xm]t  denote  a  complex  spherically 
invariant  random  vector  (SIRV)  having  zero  mean,  positive 
definite  Hermitian  covariance  matrix  R  and  characteristic 
PDF  fv(v).  The  PDF  of  x  is  given  by  [15] 

f(x)  =  n~M\R\~1h2M(q)  (1) 

where  | .  I  denotes  determinant  and 

q  —  xhr-1x 

h2M {w)  =  f0°°  v~2Mexp(- %)fv ( v)dv .  ' " 5 

Every  SIRV  admits  a  representation  of  the  form  [16]  x  = 
zV,  where  z  has  a  complex-Gaussian  PDF,  CV(0.  R),  and 
V  is  a  statistically  independent  random  variable  with  PDF 


fv(v)-  In  practice,  R  and  fv(v)  are  unknown.  For  the  pur¬ 
pose  of  this  paper,  we  assume  knowledge  of  fv  (v)  and  treat 
the  problem  of  non-homogeneity  detection  with  respect  to 
unknown  R. 

Previous  work  [1-4, 8, 12-14]  employed  the  GIP  based 
NHD  for  Gaussian  interference  scenarios.  However,  the 
GIP  based  method  relying  on  the  statistics  of  Q  =  xHR~1x 
is  unsuitable  for  SIRV  interference  scemios.  This  is  due  to 
the  fact  that  the  covariance  matrix  estimate  for  this  problem 
can  be  obtained  to  within  a  constant  of  the  covariance  ma¬ 
trix  underlying  the  Gaussian  component  of  the  SIRV.  Typ¬ 
ically,  this  constant  is  unknown  in  practice.  Consequently, 
the  PDF  of  Q,  its  moments,  and  the  threshold  setting  for  the 
goodness-of-fit  test  proposed  in  [13]  cannot  be  determined. 
Consequently,  we  seek  a  test  statistic,  which  is  invariant  to 
the  unknown  scaling. 

3.  NONHOMOGENEITY  DETECTOR  FOR 
NON-GAUSSIAN  INTERFERENCE  SCENARIOS 


where 


—  h2M^ j) 


h'2M(w)  =  =  -h2M+2(w) 


(4) 


and  qi  =  xffR-1Xj,  i  =  1, 2, ...  K .  Clearly  the  transcen¬ 
dental  nature  of  the  estimate  precludes  obtaining  a  closed 
form  solution.  Consequently,  [19]  used  the  EM  algorithm 
to  obtain  an  iterative  solution  to  the  problem.  We  adopt  the 
approach  of  [19]  for  obtaining  the  covariance  matrix  esti¬ 
mate  in  this  work.  A  derivation  of  the  covariance  matrix 
estimate  is  contained  in  Appendix  A.  It  is  shown  in  Ap¬ 
pendix  A  that  the  EM  algorithm  at  convergence  produces  an 
estimate  which  is  to  within  a  multiplicative  constant  of  the 
covariance  matrix  estimate  of  the  Gaussian  component  un¬ 
derlying  the  SIRV.  Details  pertaining  to  the  initial  start  and 
convergence  properties  of  the  EM  algorithm  can  be  found 
in  [19].  The  next  step  is  to  use  this  estimate  in  a  maximally 
invariant  decision  statistic  for  non-homogeneity  detection. 


Let  x  ~  5//2V'[0,  Rr,  fv{v)}  denote  the  complex  SIRV 
test  data  vector,  where  Ry  is  unknown.  Further,  x*,  i  = 
1, 2, ...  K  denote  iid  complex  SIRV[ 0,  Rx,  fv(v)]  target 
free  training  data.  For  homogeneous  training  data,  Ry  = 
Rx  =  R.  The  first  step  in  deriving  the  NHD  detector  for 
SIRVs  involves  obtaining  the  maximum  likelihood  estimate 
of  the  underlying  covariance  matrix.  This  estimate  is  then 
used  in  a  test  statistic  which  exhibits  maximal  invariance 
with  respect  to  the  unknown  scaling  of  the  estimated  co- 
variance  matrix.  The  resulting  test  statistic  takes  the  form 
of  a  normalized  adaptive  matched  filter  (NAMF),  which  has 
been  extensively  analyzed  in  [17, 18]  and  references  therein. 

3.1.  Covariance  Matrix  Estimation 

The  unknown  covariance  matrix  is  estimated  from  target 
free  training  data  consisting  of  independent  identically  dis¬ 
tributed  SIRVs  sharing  the  covariance  matrix  of  the  noise  in 
the  test  cell.  Maximum  likelihood  (ML)  estimation  of  the 
covariance  matrix  for  SIRVs  was  first  considered  in  [19]. 
The  work  of  [19]  showed  that  covariance  matrix  estimation 
for  SIRVs  can  be  treated  in  the  framework  of  a  complete- 
incomplete  data  problem  and  pointed  out  that  the  maximum 
likelihood  estimate  of  the  covariance  matrix  is  a  weighted 
sample  matrix.  Since  the  problem  does  not  permit  a  closed 
form  solution,  [19,20]  uses  an  iterative  method  known  as  the 
expectation-maximization  (EM)  algorithm.  More  precisely, 
let  Xj,  i  =  1,2 ,,K  denote  independent  identically  dis¬ 
tributed  training  data  sharing  the  covariance  matrix  of  the 
test  data  vector  x.  The  work  of  [19, 20]  shows  that  the  ML 
estimate  of  the  covariance  matrix  is  given  by 

1  K 

R  =  CiXiXi  (3) 

i=  1 


3.2.  Maximally  Invariant  NHD  Test  Statistic 

The  maximal  invariant  statistic  for  different  scaling  of  test 
and  training  data  is  given  by  [17] 

.  -  (5) 

[sH  R~ 1  s]  [x-^  R- 1  x] 

where  s  =  ^==[1  1  ...1]T.  Invariance  properties  of 

the  test  statistic  of  eq  (5)  and  its  geometrical  representa¬ 
tion  have  been  studied  in  [17]  and  references  therein  for  the 
case  of  Gaussian  interference  statistics  using  a  sample  co- 
variance  matrix  estimate.  In  SIRP  interference,  however, 
each  training  data  vector  realization  is  scaled  by  a  differ¬ 
ent  realization  of  V.  Consequently,  maximal  invariance  of 
the  test  statistic  of  eq  (5)  afforded  by  the  sample  covariance 
matrix  estimate  no  longer  applies.  This  is  due  to  the  fact 
that  the  sample  covariance  matrix  is  no  longer  the  maxi¬ 
mum  likelihood  estimate  of  the  covariance  matrix  for  SIRV 
scenarios.  However,  using  an  estimated  covariance  matrix 
of  the  form  of  eq  (3)  restores  the  maximal  invariance  prop¬ 
erty  of  the  test  statistic  of  eq  (5).  This  is  due  to  the  fact  that 
the  resultant  covariance  matrix  estimate  is  to  within  a  mul¬ 
tiplicative  constant  of  the  covariance  matrix  corresponding 
to  the  Gaussian  component  of  the  SIRV. 


3.3.  PDF  and  Moments  of  the  Non-Gaussian  NHD  Test 
Statistic 


The  PDF  and  moments  of  the  NHD  test  statistic  are  readily 
determined  in  terms  of  the  corresponding  quantities  of  an 
equivalent  random  variable  defined  by 


A-namf 

1  -  A-NAMF 


(6) 


NHD  Test  Statistic  PDF 


It  has  been  shown  in  [17,21,22]  that  Aeq  admits  a  represen¬ 
tation  as  the  ratio  of  an  F-distributed  random  variable  and  a 
beta-distributed  loss  factor.  In  this  effort,  we  are  interested 
in  the  PDF  of  A namf  under  the  condition  where  no  tar¬ 
get  is  present  in  the  test  data  vector  x.  Specifically,  it  can  be 
shown  from  the  work  of  [17,21,22]  that  the  PDF  of  A  namf 
is  given  by 


/a 


NAMF 


L(1  ~7)/r(7)(l  ~r)L  ^7 
[1  -  7 r]L+l 


(7) 


where  L  =  K  -  M  +  1  and  T  is  the  loss  factor  random 
variable,  whose  PDF  is  given  by 


/r(7)  = 


/3(L  +  1,M  —  1) 


7i(l-7) 


M- 2 


(8) 


The  mean  of  A  a -amf  is  somewhat  difficult  to  calculate. 
Consequently,  we  work  with  the  mean  of  Aeq  given  by 


E(  Aeg)  = 


K 


(K-M)(M-  2)' 


(9) 


The  statistical  equivalence  of  Aeq  to  within  a  scalar  of  the 
ratio  of  a  F-distributed  random  variable  and  a  beta-distributed 
loss  factor  in  that  it  permits  rapid  calculation  of  the  mo¬ 
ments  of  Aeq.  More  importantly,  it  is  extremely  useful  in 
Monte-Carlo  studies  involving  simulation  of  A  namf-  For 
homogeneous  training  data,  the  use  of  (6)  circumvents  the 
need  to  explicitly  generate  the  test  data  vector  x  and  the 
training  data  vectors  used  for  covariance  estimation.  For 
large  M  and  perforce  K.  significant  computational  savings 
can  be  realized  from  the  method  of  (6).  It  is  instructive  to 
note  that  the  PDF  of  Anamf  as  well  as  its  mean  depend 
only  on  M  and  K ,  which  are  under  the  control  of  a  sys¬ 
tem  designer,  and  not  on  nuisance  parameters  such  as  the 
true  covariance  matrix  underlying  the  interference  scenario. 
Furthermore,  for  K  — >  oc  the  mean  of  eq  (9)  converges 


to  E{Aeq) 


— — — ,  corresponding  to  the  mean  of  an 


F-distributed  random  variable. 


3.4.  Goodness-of-Fit  Tests 

Since  the  PDF  and  mean  of  An  amf  are  known,  a  formal 
goodness  of  fit  test  can  be  used  for  non-homogeneity  detec¬ 
tion  in  non-Gaussian  interference  scenarios.  In  particular, 
we  form  empirical  realizations  of  A  namf  from  each  train¬ 
ing  data  realization  using  a  moving  window  approach.  In 
this  approach  each  training  data  vector  is  treated  as  a  test 
cell  data  vector,  whose  covariance  matrix  is  estimated  from 
neighboring  cell  data  according  to  eq  (3).  We  then  test  for 
statistical  consistency  of  these  realizations  of  A  namf  with 
the  PDF  of  eq  (7).  For  this  purpose  a  convenient  type-I  error 
(typically  between  0.01  and  0.1)  given  by 

p-=JU  4~wd->  ™ 


Fig.  1.  NHD  Test  Statistic  PDF 


Fig.  2.  Type-I  Error  vs  Threshold  for  NHD  Test  Statistic 


is  chosen.  The  threshold,  77,  is  determined  by  a  numerical 
inversion  of  eq  (10).  Realizations  of  Aa tamf,  which  exceed 
77  correspond  to  nonhomogeneous  training  data.  A  second 
test  for  training  data  nonhomogeneity  is  based  on  compar¬ 
ing  each  realization  of  Aeq  with  its  theoretically  predicted 
mean  given  by  eq  (9)  and  retaining  those  realizations  which 
exhibit  the  least  deviation.  Performance  analysis  of  these 
NHD  methods  is  presented  in  the  next  section. 


4.  PERFORMANCE  ANALYSIS 

Performance  of  the  goodness-of-fit  test  with  simulated  and 
measured  data  is  presented  here.  Figure  1  shows  the  plot 
of  the  PDF  of  A  a tamf  with  K  as  a  parameter.  Observe 
that  the  variance  of  Ax  amf  decreases  with  increasing  K. 
Figure  2  shows  a  plot  of  the  Type-I  error  versus  the  thresh¬ 
old,  77,  with  K  as  a  parameter.  For  a  given  type-I  error,  the 
threshold  decreases  with  increasing  K,  in  conformance  with 
the  results  of  Figure  1 .  Figure  3  shows  the  performance 
of  the  goodness-of-fit  test  for  simulated  homogeneous  data 


Fig.  3.  Type-I  Error  vs  Threshold  for  NHD  Test  Statistic 


Fig.  4.  Type-I  Error  vs  Threshold  for  NHD  Test  Statistic 

from  the  K-distribution  [15]  with  shape  parameter  0.5  us¬ 
ing  the  covariance  estimate  of  eq  (3).  The  results  reveal  the 
lack  of  nonhomogeneity  in  that  no  realization  of  Anamf 
exceeds  the  threshold.  Figure  4  shows  the  performance  of 
the  goodness-of-fit  test  in  non-homogeneous  K-dsitributed 
clutter  with  shape  parameter  0.5.  The  K-distributed  ampli¬ 
tude  PDF  is  given  by 

aa+lra 

Jr (r)  =  2Q-ir (a)Ka~1^  r>0,/3,a>0  (11) 

where  /?  and  a  are  the  distribution  scale  and  shape  param¬ 
eters,  respectively,  K„(.)  is  the  modified  Bessel  function  of 
the  second  kind  of  order  v  and  T(.)  is  the  Eulero-Gamma 
function.  Small  values  of  a  result  in  heavy-tails  for  the  PDF 
of  (1 1).  The  corresponding  fv(v)  and  are  given  by 

rfij  (&v)2a-1exp{-/3'2v2)  0<v<oo 
h2M{w )  =  2/V^) 

Nonhomogeneity  of  the  data  is  evident  in  those  range  bins 
where  Anamf  exceeds  77.  Figure  5  shows  the  results  of  the 
goodness  of  fit  test  for  the  MCARM  data  [23]  using  acqui¬ 
sition  220  on  Flight  5,  cycle  e  for  8  channels  and  16  pulses. 


Test  Statistic  vs  Range 
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Fig.  5.  NHD  Test  Statistic  vs  Range 
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Fig.  6.  NHD  Test  Statistic  vs  Range 


The  NHD  test  statistic,  Anamf ,  and  the  threshold,  77,  are 
plotted  as  a  function  of  range.  Statistical  analysis  of  the  data 
indicates  that  the  data  is  well  approximated  by  the  Gaussian 
distribution.  This  fact  considerably  simplifies  the  analysis 
in  that  the  covariance  matrix  estimate  is  simply  the  sample 
covariance  matrix.  Non-homogeneity  of  the  training  data  is 
evident  in  those  bins  for  which  Anhd  exceeds  77.  Figure 
6  shows  the  results  from  the  selection  procedure  based  on 
comparing  the  empirically  formed  Aeq  with  its  theoretically 
predicted  mean  given  by  eq  (9).  The  data  set  used  for  this 
example  is  identical  to  the  data  used  in  Figure  5.  Relevant 
test  parameters  are  reported  in  the  plot.  We  observe  a  signif¬ 
icant  increase  in  the  number  of  deviations  from  the  theoreti¬ 
cally  predicted  mean  given  by  eq  (9).  This  is  due  to  the  fact 
that  we  are  dealing  with  a  limited  number  of  realizations  of 
the  NHD  test  statistic. 

5.  CONCLUSION 

This  paper  provides  a  rigorous  statistical  characterization 
of  the  NHD  for  non-Gaussian  interference  scenarios  which 


can  be  modeled  as  a  spherically  invariant  random  process. 
It  is  noted  that  the  NHD  statistic  admits  a  simple  represen¬ 
tation  in  terms  of  a  ratio  of  an  F  distributed  random  variable 
and  a  beta  distributed  loss  factor.  A  formal  goodness-of-fit 
test  based  on  this  representation,  which  follows  a  random¬ 
ized  F-distribution,  is  derived.  Performance  analysis  of  the 
method  is  considered  in  some  detail  using  measured  data 
from  the  MCARM  program.  The  illustrative  examples  val¬ 
idate  the  approach  taken  and  confirm  the  results.  Future 
work  would  include  extensive  performance  analysis  using 
simulated  and  measured  data  showing  the  resulting  impact 
on  STAP  performance.  The  performance  of  several  STAP 
algorithms  in  Gaussian  and  non-Gaussian  interference  sce¬ 
narios  has  been  considered  in  [18].  Future  work  will  ad¬ 
dress  performance  of  the  methods  treated  in  [18]  with  suit¬ 
able  NHD  pre-processing. 
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is  either  z u  Vi,i  —  1,2, . . .  ,K  or  x*,  Vi,  i  =  1,2,...,  iV. 
However,  the  observed  data  Xj,  i  =  1,2, ....  K,  contains  no 
explicit  information  about  Vi,  i  =  1, 2, . . . ,  K  and  thus  con¬ 
stitutes  the  incomplete  data.  The  complete  data  likelihood 
function  is  given  by  the  joint  PDF  of  Xj,  Vj,  i  =  1, 2, . . . ,  K, 
which  is  expressed  as 

K  K 

5c[x,  Vi\R] = n /tew  n  /w  <i4> 

i=l  i=l 

Taking  the  natural  logarithm  of  (14)  yields  the  complete- 
data  log-likelihood  function  of  the  form 

K 

L[X,Vi\R]  =  —KMlog(it)  -  Klog{ |R|)  - 

K  *=1 

+Y^lo9[v~2M  }{vi)]. 

i=  1 

„  (15) 

Note  that  given  an  initial  estimate  of  R  denoted  by  R,  the 
quantity, 

E{log[v~2Mf(Vi)]\R}  (16) 

depends  only  on  R  and  not  on  R.  Consequently,  the  rele¬ 
vant  terms  for  the  maximization  over  R  are  given  by 

K 

L:[X,  Vi\R)  =  — iTZop(|R|)  -  ^vf2  (17) 

»=i 


[24]  J.  P.  Burg,  D.  G.  Luenberger,  and  D.  L.  Wenger,  “Es¬ 
timation  of  structured  covariance  matrices,”  Proceed¬ 
ings  of  the  IEEE,  vol.  70,  No.9,  pp.  963-974, 1982. 

Appendix  A:  EM  Algorithm  for  Covariance  Ma¬ 
trix  Estimation 

We  discuss  the  maximum  likelihood  estimation  of  the  SIRV 
covariance  matrix  in  this  appendix.  Let  X  denote  a  data 
matrix,  whose  columns  Xj  i  =  1,2,...,  iV  are  indepen¬ 
dent  identically  distributed  target-free  training  data  vectors, 
which  are  distributed  as  SIRV[0,  R^.,  /y  (v)].  The  likeli¬ 
hood  function  for  estimating  R  if  given  by 

K 

<z[X|R]  =  n  tt-m|R  r1  W«i).  (13) 

i=l 

Direct  maximization  of  the  likelihood  function  of  (13)  over 
R  is  rendered  difficult  due  to  the  fact  that  there  is  miss¬ 
ing  information.  Consequently,  it  is  helpful  to  treat  the 
problem  in  the  context  of  a  complete-incomplete  data  prob¬ 
lem  [19],  Recall  from  the  representation  theorem  for  SIRVs 
[16]  that  Xj  =  ZjVj,  where  z*,  i  —  1,2,...,  AT  are  sta¬ 
tistically  independent  CN( 0,  R)  random  vectors,  and  V 
i  =  1, 2, . . . ,  K  are  statistically  independent  random  vari¬ 
ables  with  PDF  fv(v).  For  this  problem,  the  complete  data 


The  missing  data,  Vi,  i  =  1, 2, ...  K,  are  assumed  to  be 
missing  at  random  (MAR)  [19].  Consequently,  given  an  ini¬ 
tial  estimate  of  R  denoted  by  R,  the  complete  data  sufficient 
statistic  [19]  is  given  by 

c*  =  E[V~2 |R,  Xi].  (18) 

Note  that  f(v i)  =  f(v 2)  =  . . .  =  f(vK)  =  fv(v)  (since 
Vi,  i  =  1, 2, ...  K  are  independent  identically  distributed 
random  variables).  Therefore, 


f(xi\vi,  R )fv(vj) 
/xi|R/Xi  W 


However, 


(19) 


f(xi\vi,  R)fv(vj)  _  Vj  2Mexp(qivl  2)f{vj) 
/xi|R.(X*l^  ^2  M(Qi) 

Consequently, 


<k  =  E[Vr2  |R,  Xi]  =  - 


^2  m(^) 
h2M(qi ) 


(20) 

(21) 


Having  specified  the  complete  data  sufficient  statistic,  we 
seek  the  maximization  of  (17).  For  this  purpose,  we  re¬ 
produce  the  following  matrix  differentiation  identities  from 
[24]. 

S[  R-1]  =  -R-1<5[R]R-1 
<5[/op|R_1|]  =  -tr{R-15[R]}  K  } 


Further,  we  recognize  that  qi  =  X^i  Con¬ 

sequently, 


K 

<5Li[X,  Vf|R]  =  KtriR-^lRjj-trlR-'SlRlR-'Y,^** 

i=  1 

(23) 

Maximization  of  (17)  results  from  setting  (23)  equal  to  zero. 
Therefore,  the  maximum  likelihood  estimate  of  R  is  given 
by 

1  K 

R=— ^CiXjxf.  (24) 

i=l 

Clearly  the  transcendental  nature  of  the  estimate  precludes 
obtaining  a  closed  form  solution. 

In  summary,  the  EM  algorithm  for  the  problem  of  co- 
variance  matrix  estimation  considered  here  consists  of  the 
following  steps. 

1.  E-Step:  Given  an  initial  estimate  of  R  denoted  by  R, 
calculate  cz  for  i  =  1, 2, . . . ,  K. 

2.  M-Step:  Calculate  R  =  ^  i  c{x.ixf . 

3.  Iterate  until  convergence.  Convergence  is  determined 
through  a  suitable  error  criterion. 

In  this  paper,  the  convergence  criterion  is  an  error  of  10“6 
defined  to  be  the  absolute  value  of  the  difference  between 
the  values  of  R  resulting  from  two  successive  iterations. 
Convergence  of  the  algorithm  is  dictated  by  the  choice  of 
the  initial  estimate  of  R.  Any  positive  definite  Hermitian 
matrix  is  suitable  for  the  initial  estimate  of  R.  Two  choices, 
that  arise  readily  are  the  M  x  M  identity  matrix  and  the 
sample  covariance  matrix  given  by  S  =  £i=1  xixf/-  We 

employ  the  latter  choice  in  this  paper  due  to  the  fact  that  it 
yields  faster  convergence. 

The  simulated  data  examples  considered  in  this  paper 
involve  the  calculation  of  the  modified  Bessel  function  of 
the  second  kind  for  specifying  /i2m(-)  and  its  derivative. 
Numerical  errors  in  their  calculation  for  a  =  0.1  tend  to  be 
rather  large.  Consequently,  convergence  of  the  algorithm  is 
extremely  slow  for  a  =  0.1. 


